© The Authors 2024

1. Introduction

Type Ia supernovae (SNe Ia) have been used for a long time as standardisable candles in the extraction of cosmological parameters (Perlmutter et al. 1999; Riess et al. 1998), and the largest correction in the standardisation scheme is related to the extinction by the host galaxy. However, a significant uncertainty remains on this reddening correction, which affects the potential of SNe Ia for the accurate determination of cosmological parameters. Here, we try to answer some remaining questions related to the use of the semi-empirical extinction formulae derived for our galaxy by Fitzpatrick (1999). A widely employed standardisation method using SNe Ia is provided by SALT2 (Guy et al. 2007; Betoule et al. 2014). This latter incorporates purely empirical correlation between colour and magnitude with the general form m_B = M₀ + αx₁ + βc, where x₁ characterises the shape of the light curve, and c ∼ (m_B − m_V) the colour of the supernova (SN). The parameters α and β are tuned so as to improve the standardisation of the SNe Ia in a given sample. This description is simple, but explicitly abstains from separating the intrinsic and reddening contributions to the observed colour. The reddening is caused by a variable mix of gaseous contributions (H₂, H, CH, etc.) and dust grains of variable sizes. Such a complex mixture would not be expected a priori to lead to a simple extinction law, but detailed work by several authors extending over 50 years (Rieke & Lebofsky 1985; Cardelli et al. 1989, hereafter CCM; Fitzpatrick 1999, hereafter FM99) led to universal formulae –depending on a single parameter R_V– related to the cross sections of the diffusion centres (molecules or dust grains) for photons. R_V is operationally defined as the ratio of the extinction in the V bandpass A_V to the reddening indicator E′∼E(B − V), either measured or adjusted by a fit to the observations.

A detailed study by Schlafly et al. (2016, 2017) confirmed a directional as well as distance variability of R_V in our galaxy (from 3 to 3.75), and provided confirmation of the results of FM99. All attempts to determine the parameter R_V from SNe assume that the extinction formulae derived for our galaxy apply directly to all host galaxies, leaving the possibility of different values of R_V, as the averaging which occurs when considering a sample of host galaxies differs from that performed when deriving the extinction formulae from stars of the Milky Way. The same universal extinction formula is assumed to be valid for the host galaxies of the SNe Ia, but previous investigations actually relied on a combination of photometric measurements with slit spectroscopy, which lacked the spectrophotometric information provided by the Nearby Supernova Factory (SNfactory, Aldering et al. 2020) collaboration. The contribution of the intrinsic variability was not directly monitored, and the applicability of the extinction formula to other galaxies could not be checked accurately. All the groups restrict their extinction analysis to redshifts larger than 0.01 in order to avoid a significant contribution of peculiar velocities.

2. Previous investigations

There have been many attempts to separate the intrinsic and extinction colour components of SNe, and we describe a few of them here. In many instances, an a priori distribution of the extinction is assumed; for instance an exponential distribution with a cutoff. This introduces a bias in the determination of the extinction of each SN, as the reddening involves an average over many different geometrical configurations of the SN with respect to its host. The exponential behaviour is a simple but unrealistic assumption. In the present work, an extinction colour and an intrinsic colour are extracted for each SN. The intrinsic colour contribution to E(B − V) is frequently derived from the light-curve shape, whether it is the SALT2 x₁ or the stretch parameter. In our case, an extra intrinsic component is introduced (for each SN).

Lira et al. (1998) and later Phillips et al. (1999) noticed that the colour evolution of all SNe Ia in the 30–90 days post V-maximum is universal, which allowed them –by selecting SNe in E or S0 galaxies (without dust)– to derive a universal intrinsic SN Ia colour E(B − V)₀ at any reference date in the interval from 30 to 90 days; the magnitude error quoted is 0.05. The galaxy reddening E(B − V)_tail at late epochs is then found by subtracting the intrinsic colour from the observed colour. After selecting a sample of SNe with low extinction (E(B − V)_tail < 0.06), a time-dependent correction provided by the light curve allows them to evaluate the intrinsic value E(B − V)_max, 0 at maximum B luminosity.

Another photometric technique was used by Wang et al. (2003), who analysed the colour–magnitude diagrams using the data from Hamuy et al. (1996) and Riess et al. (1998). The B magnitude varies linearly as a function of colour in the 10–30 day period post B maximum. The scatter of the residuals from the straight line is typically 0.05. It is claimed that extinction does not affect the shape of these diagrams (an approximation, as the spectra of the SNe evolve). The magnitude $M_{BV}^{\prime }$ of the SN Ia at colour B − V = 0.6 is used as a reference value for the standardisation. To obtain a ‘reference’ magnitude at maximum, the linear behaviour is extended to B − V = 0, although some SNe Ia show a ‘bump’ departing from linearity. In the absence of a such a bump, the comparison with the observed $M_{\max }^B$ allows a value of the reddening E(B − V) to be derived, which has a dispersion of about 0.1 mag with respect to the determination of Phillips et al. (1999) for the same SN Ia.

Jha et al. (2007) also relied on photometric bandpasses and K corrections, but took into account the colour evolution along the light curve of the SNe using the MLCS2k2 software (Multicolor Light Curve Shapes). A preliminary fit extracts the intrinsic Gaussian fluctuation (with σ = 0.049, and an average value of zero) and the exponential reddening distribution at date +35 days (with a colour decay constant of τ_E = 0.138). In each UBVRI bandpass, the same reddening distribution is used in the modelling of the light curve. Each of the five light curves is described by five parameters: date of peak B luminosity, distance modulus, time evolution, a unique reddening scale (A_V) (with the B  −  V shape found in the previous step), and R_V. A prior R_V = 3.1 with σ = 0.4 is assumed for the CCM law (the parameter R_V is needed to convert A_V into an extinction in different filters). These latter authors do not claim to describe the intrinsic distribution at maximum luminosity, nor to measure R_V.

The multi-colour light curves of a sample of 80 SNe measured by different groups are analysed by Nobili & Goobar (2008). The light curve shape is characterised by its stretch, that is, the time dilatation factor s(sn) bringing the mean light-curve width to the SN light-curve width. The stretch is similar to the x₁ variable mentioned earlier for SALT2. Given two bandpasses, X and Y, and stretch s: X − Y = b_XY(t)+a_XY(s − 1)+c_XYE(B − V). The data are K-corrected for the changes of rest-frame bandpasses with redshifts, and E(B − V) is the average over all epochs of the B − V colour curve of each SN. R_V is first extracted from the measured coefficients c_XY found (e.g. c_BV = R_B − R_V). The best fit is found (with the CCM law) for R_V = 1.01 ± 0.25. When the excess dispersion of the residuals along the light curve after correction for the extinction is assigned to the intrinsic colour contribution, the same analysis finds R_V = 1.75 ± 0.27 using CCM. This difference strongly suggests that the intrinsic colour components should be taken into account.

The ‘twin’ SNe SN2014J and (unextincted) SN2011fe were compared in the optical and near-infrared (NIR) range by Amanullah et al. (2014) from data collected by the HST (UV bands), at Mauna Kea (NIR), and NOT (optical). The comparison of the 12 light curves was performed using the extinction formula of FM99. Amanullah et al. (2014) find E(B − V) = 1.37 ± 0.03, R_V = 1.4 ± 0.1.

The Carnegie Supernova Project (Burns et al. 2014) emphasises the use of a ‘pseudo-stretch’ s_BV, which describes the time dependence of the B − V colour instead of the usual ‘stretch’, which relates directly to the light curve. This pseudo-stretch allows a prediction of the maximum B − V colour. The analysis proceeds to analyse the eight observed colours B − m_i in bandpasses (u, g, r, i, V, Y, J, H) of each SN as a sum c_i(sn) = P(s_BV)+ΔA_i(E(B  −  V),R_V), where P(s_BV) is a second-order polynomial, and ΔA_i is the predicted reddening for colour c_i from the extinction formula. From a set of 75 SNe, Burns et al. (2014) find R_V = 2.14 ± 0.16. With one intrinsic colour and one (measured) host extinction per SN, and without any assumption on the distribution of the extinction, their method and their result are relatively similar to ours. In Sect. 9, we suggest using the SALT2 variable x₁ in the way s_BV is used by these authors.

Amanullah et al. (2015) include photometric measurements of seven SNe Ia between 0.2 and 2 microns from date −10 to +50. A mean spectral template is used at each date to evaluate the extinction correction in all filters with a standard CCM extinction law where E(B − V) and R_V are parameters to be adjusted to the observed colour. Amanullah et al. (2015) assume an intrinsic U − V dispersion of 0.1 mag (compatible with our own spectrophotometric measurements), and a larger dispersion of 0.3 mag for wide-band filters at lower wavelengths. The different intrinsic colours have a mean value of zero and are supposed to be uncorrelated at any given date, but fully correlated at different phases. The light curves are then fitted in all filters as in Amanullah et al. (2014) to extract R_V and E(B − V) for each SN. The values of R_V for the different SNe vary from 1.4 ± 0.1 to 3.8 ± 1.5.

Mandel et al. (2017) analysed a sample of 250 nearby SNe Ia (0.01 < z < 0.10). They adapted the SALT2 formalism (Betoule et al. 2014) at maximum B light to include an additional intrinsic contribution of the intrinsic colour c^int to the B magnitude as well as an explicit reddening R_BE_s. The intrinsic B band absolute magnitude is given by $M^{\text{int}}=M_0^{\text{int}}+\alpha x^{\text{int}}+{\beta }_{\text{int}}{c}^{\text{int}}+ϵ$ with a contribution of the intrinsic colour c^int to the B magnitude in addition to the reddening R_BE_s. The intrinsic light-curve shape parameter x, similar to SALT2 x₁, is assumed to have a normal distribution x^int ≈ 𝒩(−0.40, σ = 1.2), the intrinsic colour is $c^{\text{int}}=c_0^{\text{int}}+{\alpha }_c^{\text{int}}{x}^{\text{int}}$ with $c_0^{\text{int}}\approx \mathcal{N}(-0.06,\sigma =0.06)$, and the extra (grey) magnitude dispersion ϵ ≈ 𝒩(0.,σ = 0.1) is introduced. The ‘grey’ fluctuation is the dispersion in the magnitude of SNe that have identical or quasi-identical spectral shapes after reddening correction; it is bound to be achromatic by definition, and its origin is not understood. The intrinsic parameters of SN s are thus $(M_s^{\text{int}},c_s^{\text{int}},x^{\text{int}})$. The host galaxy reddening E_s is positive, with an exponential dependence on colour exp(−E_s/τ), as in Jha et al. (2007), with a magnitude scale defined by R_B. The model is adjusted to the outcome of a SALT2 fit (m_B, c, x₁); the Hubble residuals are 0.16 mag. The 11 parameters of the model are obtained by a hierarchical Bayesian method. The observed values, which include reddening and distance modulus, are $(m_s,c_s^o,x_s)$. The main results of the fit in Mandel et al. (2017) are: $c_0^{\text{int}}=-0.061\pm 0.012$, σ^int = 0.100 (intrinsic magnitude dispersion), ${\sigma }_C^{\text{int}}=0.065\pm 0.008$ (intrinsic colour dispersion), R_B = 3.73 ± 0.31, τ_E = 0.069 (extinction decay rate). The dependence of the magnitude with intrinsic colour is observed with β^int = 2.25 ± 0.25. The dispersion of the intrinsic B − V colour found by these latter authors is four times larger than our result.

Thorp et al. (2021) used photometric data to investigate the potential dependence of the reddening as a function of host mass. Their model includes the Gaussian grey magnitude fluctuation as well as a spectral variability function (largely spectral lines) with an amplitude scale as a parameter. The R_V values of individual SNe are drawn from a truncated Gaussian distribution (R_V > 0.5). The Hubble residuals are improved from 0.16 (SALT2) to 0.12. Thorp et al. (2021) finds no significant difference in R_V between the large galactic mass and low-mass samples. The full sample averages to R_V = 2.70 ± 0.25.

As in Jha et al. (2007), Brout & Scolnic (2021) describe the observed colour as a sum of intrinsic and extinction colours, with a Gaussian distribution of the intrinsic colour, and an exponential shape of the dust. Intrinsic and extinction colour components are added into the observed colour. Their model also includes a Gaussian R_V distribution. When splitting their sample into low-mass and high-mass galactic hosts, Brout & Scolnic (2021) find two different values of R_V, with R_V = 2.0 ± 0.25 (low mass) and R_V = 3.0 ± 0.4 (high mass).

A recent study by Wojtak et al. (2023) also investigates the impact of the host mass on the standardisation of SNe Ia within a Bayesian approach using multi-filter light curves of SNe Ia. These authors introduce a shape parameter for the extinction distribution rather than an exponential with a cutoff and find strong evidence for two populations of SNe. The two populations, with probabilities 0.38 and 0.62, differ in their intrinsic colour and their reddening distribution. Wojtak et al. (2023) introduce an (unaccounted) intrinsic magnitude dispersion of 0.11 mag. The authors also define two intrinsic colour dispersions for the two populations: ${\sigma }_{\text{int}}^1=0.076\pm 0.011$ and ${\sigma }_{\text{int}}^2=0.043\pm 0.013$. Both these values are larger than our own determination, which includes more spectroscopic information.

The only spectrophotometric data available are provided by the SNfactory with the SNIFS instrument. A first attempt by Chotard et al. (2011, hereafter C11), analysed the magnitudes of 76 SNe Ia at B maximal luminosity in the five synthetic top-hat filters UBVRI. The magnitudes were corrected for the Si II and Ca II H&K equivalent widths, and the extinction correction was then derived from the CCM law with scale A_V. A value of R_V = 2.8 ± 0.30 was obtained. The extinction scale A_V was a free parameter adjusted to the data for each SN in this latter study. The magnitude uncertainties were estimated from the measurement errors (dominated by the flux calibration uncertainty) without including the (unaccounted) grey fluctuation of 0.11 mag. Colours are better measured than magnitudes as some systematic uncertainties cancel out, such as the absolute flux calibration uncertainty of 0.03 mag, and (largely) the error on the redshift (measurement and peculiar velocity). The present work is an extension of this latter investigation, but with an increased sample and several improvements.

Huang et al. (2017) uses the twin SNe SN2012cu (with a significant reddening) and SN2011fe (without extinction) measured by the SNfactory to derive the parameter R_V by comparing their magnitudes at different epochs. Although the spectra are very similar after extinction correction (the two SNe are ‘twins’), some spectral bins near absorption lines differ. These bins are averaged by a Gaussian convolution, and deweighted so that the χ² in these regions is unity (a ‘floor’ error of 0.03 is added); the residual fluctuations of the deweighted residuals are of the order of 0.01 mag, and the value of R_V is R_V = 2.952 ± 0.081 for SN2012cu. Léget et al. (2020) and Boone et al. (2021) used the SNfactory spectrophotometric data to show that the largest spectral features occur with a remarkable intrinsic variability in particular around Ca II H&K and Si II.

While Léget et al. (2020) singles out the correlation of spectral lines with magnitudes, Boone et al. (2021) introduced three intrinsic variables ξ_i, which provide an accurate description of the full spectra for 173 SNe Ia (after selection). The formula of FM99 is used to correct for reddening. As in Huang et al. (2017), the parameter R_V is found by deweighting the spectral ranges with large intrinsic variability by a factor 1/σ(λ) at each wavelength λ, with the result R_V = 2.40 ± 0.16 (the method is called Read Between The Lines). The amplitude of the residual intrinsic variability of magnitudes away from spectral lines after reddening correction is generally of the order of 0.08 mag, but it is as low as 0.02 mag between 6600 and 7200 Å.

A large uncertainty remains today on the (mean) value of the parameter R_V, which plays a substantial role in the standardisation of SNe. Part of the discrepancy arises from underestimated systematic measurement errors on the path from detector to magnitude: the absolute calibration, the bandpass wavelength dependence, the redshift accuracy, the peculiar velocities, and the corrections for galactic reddening (host and ours). Another source of variability lies in different definitions of R_V as a parameter of an extinction formula, or as A_V/E(B − V), where A_V and B − V depend on the spectral distribution of the stars used (HD stars and SNe) and on the rest-frame bandpasses selected. Peculiar velocities are part of the redshift accuracy and influence the rest frame bandpasses as well as the absolute magnitude.

The spectrophotometric quality of the SNfactory data (Aldering et al. 2020) helps alleviate some of these uncertainties, and the difference in R_V between Huang et al. (2017) and Boone et al. (2021) –who average over 173 SNe– as well as the results of Amanullah et al. (2015) –who use a somewhat less homogeneous sample– may well point to an actual physical distribution of the value of R_V. Here, to quantify the smooth spectral evolution of the reddening, we favour the use of colours rather than magnitudes, as the latter exhibit an ‘unaccounted’ grey fluctuation of 0.08–0.10 mag from one SN to another, and the evaluation of the ratio of the correlated to uncorrelated colour variabilities, including the intrinsic contribution, is not straightforward.

We assume in the present work that the reddening of SNe by the host galaxy can be factored out and consider the observed colour fluctuations (from one SN to another) as a sum of intrinsic and extinction contributions, as done previously by Jha et al. (2007), Mandel et al. (2017) or Brout & Scolnic (2021). As Si IIλ6355 Å and Ca II H&Kλ3945 Å are dominant components of the spectral intrinsic variability, we first consider their contribution –as in C11–, but substituting Si IIλ4131 Å for Si IIλ6355 Å as the dependence of the B magnitude on the equivalent width of Si IIλ6355 Å is quadratic. We shall assume in addition that there is a single dominant source of intrinsic variability on top of the Ca II H&K and Si II absorption lines, so that the intrinsic variation of the different ‘colours’ is linked by three ‘intrinsic’ couplings. For each SN, we can then extract the extinction and intrinsic colour components using the extinction formula as leverage, and we revisit the determination of R_V using the correlations of these colour variations.

There are three advantages to using the particular data sample chosen here: the high-quality spectrophotometric measurement, with well-measured spectral lines and subpercent errors on colours, and the sizeable homogeneous sample of 165 SNe, which allows us to use small colour fluctuations with respect to a mean template rather than the absolute values of these colours.

3. Sample of supernovæ

The data used were gathered by the SNfactory collaboration using their SuperNova Integral Field Spectrograph (SNIFS, Lantz et al. 2004), an automated instrument optimised for the observation of point sources on a structured background, with a spectral resolution of 0.25–0.30 nm and a good efficiency from 330 to 900 nm. SNIFS consists of a multi-filter photometric channel used to monitor the transmission in non-photometric nights and provide an image for guiding, a lenslet integral field spectrograph covering a field of view of $6\stackrel{″}{.}4\times 6\stackrel{″}{.}4$ with a grid of 15 × 15 spaxels, and an internal calibration unit (continuum and arc lamps). A more complete description of SNIFS, its operation, and the data processing can be found in Aldering et al. (2006); updated in Scalzo et al. (2010). The measured wavelength range in the SN rest-frame extends from 3300 to 8400 Å. It is divided into five logarithmically distributed top-hat filters UBVRI, with central wavelengths at 363.9, 438.6, 528.7, 637.4, and 768.4 nm. The light curves are reconstructed in these five synthetic rest-frame bandpasses, and the present study is limited to the spectrum closest in time (within a window of ±2.5 days) to the maximal B luminosity, as found from a SALT2 fit (Betoule et al. 2014) to the light curves. The redshift range covered by our measurements extends from z = 0.02 to z = 0.11. Although there are only four independent colours, we consider all pairs of filters in the present analysis, which leads to ten possible filter combinations. The initial sample consists of 172 SNe Ia, which were selected as those yielding an acceptable SALT2 fit (Betoule et al. 2014), allowing us to define the date of maximum light: at least 5 nights of observations, nMAD¹ of residuals < 0.12 mag, and a satisfactory phase coverage: at least four epochs from −10 to +35 days from maximum light, with one epoch between −10 and +7 days, and one between +7 and +20 days. In addition, we required in this work one spectrum within 2.5 days of maximum light. We eliminated seven SNe Ia where the subtraction of the host galaxy signal in our version of the data processing left a brightness gradient in either the B channel of the spectrograph (up to a wavelength of 5000 Å) or the R channel (above 5000 Å) over the 225 micro-lenses larger than 0.05 mag. The analysis is performed on the remaining 165 SNe Ia from the SNfactory. This sample is almost the same as that described in Aldering et al. (2020). As explained in the following section, the SALT2 magnitudes for our five synthetic bandpasses are not used; only the shape of the light curve near maximum light is used, as provided by the SALT2 model after the determination of x₀ and x₁ from the data.

Magnitudes at maximum light and errors. We require the magnitudes produced by the SALT2 fit in each bandpass to be consistent with the data near peak. From each rest-frame spectrum within a ±2.5 day window around the B-band maximum light, the distance-corrected magnitudes (brought back to an arbitrary common redshift of 0.05) are derived by integrating the photon count in each of the five UBVRI top-hat bandpasses described in Table 1. The index T in Table 1 is used to distinguish the bandpasses from Bessell filters, and is dropped later in the text. The magnitude in each filter is rescaled to an interpolation of the actual observations as described in C11. The SALT2 fit extracts the date of the B maximum, the values of x₀ (absolute magnitude scale) and x₁ –which relate to the light curve–, and the colour c, which is very close to B − V at maximum B light. These parameters are then used to integrate the SALT2 spectral templates over our bandpasses, and the corresponding SALT2 light curves are reconstructed at the observation dates, for each SN, according to the values of $M_F^{0S}$, x₁ and c parameters found in the light-curve fit.

In each bandpass, the shape of the SALT light curve near maximal B luminosity is then used as an interpolating curve in each filter and scaled as described in Eqs. (1) and (2) to fit our synthetic photometry from spectra within 5 days of maximum luminosity, as in C11. The same integration over the bandpass is applied to the SALT2 spectral template so that the corresponding SALT2 light curves are reconstructed at the observation dates, for each SN, according to the values of x₀, x₁, and c parameters found in the light curve fit. The SALT magnitudes $m_{F,p}^S$ reconstructed from the light curve in each bandpass F are then shifted at all dates p according to Eq. (2) so as to match the observed distance-corrected magnitude m_F, p. The shift ϵ_F is found by minimising

$$\begin{array}{c}\hfill {\chi }_F^2={\displaystyle \sum _{p=1}^N\frac{{(m_{F,p}-m_{F,p}^S-{ϵ}_F)}^2}{{\sigma }_{F,p}^2}.}\end{array}$$(1)

The values of ϵ_F are typically about 0.01 mag within the 2.5 day window around maximum light t_max, with a spread of 0.02 to 0.05 mag around the SALT2 light curve for UBVR, and somewhat larger for I. The rescaled magnitudes m_F at maximum B luminosity are then obtained by shifting the SALT maximum $m_{F,\max }^S$ by this (averaged) ϵ_F:

$$\begin{array}{c}\hfill m_F=m_{F,max}^S+{ϵ}_F.\end{array}$$(2)

Our magnitudes are then largely independent of the SALT colour law and of the SALT light-curve parameters. Given the measurement errors of the magnitudes for the different filters at each date, and the covariance of the measurement errors between the different bandpasses (mostly a grey absolute calibration factor and a redshift uncertainty from peculiar velocities), we can then estimate the covariance matrix of the magnitudes in the different filters at maximum B luminosity. The rescaled SALT2 magnitudes at maximum are corrected for the extinction in our galaxy.

4. Methodology

The data analysis presented here starts from the processing described in Aldering et al. (2020). It is followed by an evaluation of the magnitudes at maximum luminosity in each bandpass, and a measurement of the equivalent widths of the significant spectral features as performed in C11 and Chotard (2011). The first step of the analysis, as in the earlier work of C11, is to subtract in Sect. 4.1 the intrinsic colour variability correlated to Si IIλ4131 Å and Ca II H&Kλ3945 Å signals using their equivalent widths ew^Si and ew^Ca. What is referred to as ‘colour’ in our study –except in this first step– is actually the ‘colour difference’ between the measured colour and the average value of the same colour over our sample of 165 SNe after subtraction of the contribution of these two lines. It differs from the standard colour excess, measured with respect to a sample of SNe without extinction by a constant, and can be negative. This constant offset (different for each filter) has no impact on any of our results. Each colour difference e_i is the sum of two components, an intrinsic part I_i and an extinction part X_i. In our method, X_i is also an extinction variation with respect to its average value and may assume negative values. For each filter, it also differs from the full extinction by an (irrelevant) constant. We then express the ten colour differences observed (after correction for ew^Si and ew^Ca) in terms of the following relations (without any summation of the indices):

$$\begin{array}{c}\hfill e_i(n)=I_i(n)+X_i(n)+{ϵ}_i(n),i=(1\dots 10),\end{array}$$(3)

and for each colour i:

$$\begin{array}{cc}& X_i(n)={\delta }_{\mathit{ij}}(R_V)X_j(n)\hfill \\ \hfill & I_i(n)={\gamma }_{\mathit{ij}}{I}_j(n).\hfill \end{array}$$(4)

The coefficients γ_ij and δ_ij are common to all SNe, labelled by index n; ϵ_i is the residual associated with the description of the colour difference e_i within the present model. While γ_ij relates intrinsic colours and does not depend formally on R_V, its determination will be weakly dependent on the value of R_V assumed. The previous relations do not involve convolutions with measurement errors, as these ‘colour’ errors are small enough with respect to the physical quantities I_i and X_i to be ignored at this stage.

Our notations are summarised in Table 2. The colours are numbered from zero to 9 in the sequence (U − B, U − V, U − R, U − I, B − V, B − R, B − I, V − R, V − I, R − I), and the differences in their measurement accuracies is ignored. The coefficients δ_ij(R_v) relating the extinction colour components X_i are derived from the extinction formula as described in Appendix A, while the 45γ_ij, which relate the intrinsic colour variations I_i, can be expressed as a function of three of them, selected (arbitrarily) to be γ₁₀, γ₂₀, γ₃₀. Thus, the proposed model is characterised at this stage by four parameters: R_V, and the three intrinsic coefficients γ₁₀, γ₂₀, γ₃₀. However, below, we find it necessary to introduce four corrections to the extinction formula in the bandpasses BVRI, which results in a total of eight parameters. For presentation purposes, an extra parameter (not required by our analysis, where it is unity) will be introduced for the extinction scale. For each SN, two coordinates are derived from the data using the previous parameters: X₀(n) and I₀(n). We call them coordinates rather than parameters as they are computed with a fixed algorithm from the measured spectra, without any tuning, once the parameters involved in our modelling of the full sample have been found. The extinction analysis proceeds in eight stages, with the successive application of different fits. To assess the validity of the model used, we have chosen to control each successive stage and apply judgement, rather than combine the different χ², which occur along the path into a single χ² or a likelihood:

1. We subtract the contributions of the silicon and calcium equivalent widths ew^Si and ew^Ca to the colour fluctuations to ensure they are small in Sect. 4.1, allowing a linear approximation.

2. The extinctions $A_F^{\prime }$ in filter F derived from the extinction formula 11 are rescaled to A_F (Sect. 4.2.2) to ensure mathematical consistency. The colour coefficients δ_ij of Eq. (3) are invariant in this rescaling.

3. For a given set of γ_ij, the intrinsic (I_i) and extinction (X_i) colour components are found by a least square fit to the observed colours (Sect. 5.2).

4. The values of the intrinsic colour correlation coefficients γ_ij are found by iteration, and these are tuned in Sect. 5.3 by minimising the scatter of the output ratios of the 10 × 9/2 = 45 pairs of colour-ratios found from the colour fit in Eq. (29).

5. The extinction formula (FM99) is modified so as to suppress any correlation between the residuals ϵ_i derived from the fit to the extinction colour. This is achieved in Sect. 5.4 by minimising ${\chi }_{\mathit{Dex}}^2=\sum _i{(d{ϵ}_i/d{X}_i)}^2$. The consistency rescaling is then applied to the corrected coefficients.

6. Iterations are performed over steps 4 and 5 on the values of the filter extinction corrections and the intrinsic coefficients γ_ij, until neither ${\chi }_{Dex}^2$ nor the scatter can be improved by a change of 10⁻³ of the extinction corrections or the coefficients (γ₁₀, γ₂₀, γ₃₀).

7. The previous steps do not depend on an overall scale factor applied to the (filter-colour) extinction coefficients δ_Fi. For each value of R_V, this scale factor is unity within our model, but it is introduced in almost all extinction corrections under the name of A_V, where the extinction at wavelength λ is A_Vϕ(λ), with the function ϕ standing for the extinction formula. To compare with models that do not enforce the normalisation of the extinction, we introduce in Sect. 7 a (illegitimate) scale parameter s multiplying the extinction magnitude A_F, so that the extinction in bandpass i will now be s δ_i4 X₄, where colour 4 ≡ B − V (with our top-hat filters).

8. It is then observed in Sect. 7 that there is a unique value of R_V for which the extinction-corrected magnitudes in the other filters UBRI do not depend on the extinction colour. This provides our determination of R_V.

This derivation of R_V relies on the suppression of the (well-measured) colour dependencies on R_V, rather than on the fit of magnitudes involving a χ² (subject to the ‘grey’ fluctuation of 0.08–0.12 mag and other correlated effects, such as modelling and calibration). There is no attempt at any stage in the previous sequence to minimise the colour residuals or the residuals of magnitudes in the Hubble diagram.

4.1. Calcium and silicon contributions

Nugent et al. (1995) showed that the impact of Ca II H&K and Si II absorption lines can be characterised by their equivalent widths and that they correlate strongly with magnitudes, colours, and light-curve shapes. This section addresses a mathematical issue: how to suppress the part of the colour variation that depends on ew^Si and ew^Ca in the presence of a correlation between these two widths. C11, and other earlier works referenced therein, observed that there is a strong correlation between the equivalent widths ew^Si of Si IIλ4131 Å, ew^Ca of Ca II H&Kλ3945 Å, and the different colours. The equivalent widths ew^Si and ew^Ca are measured as in Bronder et al. (2008), using an algorithm described in C11, and in more detail in Chotard (2011): following Nugent et al. (1995), a spectral range centred on the spectral line is defined from λ₁ to λ₂, and a reference level f_c(λ) is defined by a straight line connecting the two spectral values associated to λ₁ and λ₂. The actual spectrum observed is s(λ). The equivalent width is

$$\begin{array}{c}\hfill ew={\displaystyle {\int }_{{\lambda }_1}^{{\lambda }_2}(1-\frac{s(\lambda )}{f_c(\lambda )})d\lambda .}\end{array}$$(5)

The errors on the equivalent widths are derived using a simulation that takes into account the photon statistics as well as the algorithm used to select the boundaries in the integration over the line width. As these two Ca and Si spectral features have a strong colour contribution over the whole spectral range, it is advantageous to treat them separately from the extra-intrinsic colour and take advantage of the detailed spectral information. The analysis described in the following sections refers to the intrinsic colour component remaining after subtraction of the effect of Si IIλ4131 Å and Ca II H&Kλ3945 Å on the ten different colours². The small (formal) variations $\delta c_i^{\prime }$ of the ‘observed’ rest-frame colour $c_i^{\prime }$ around its mean value are assumed to depend linearly on the variations (δew^Si,δew^Ca) of the equivalent widths (with other contributions treated as noise).

The subtraction of the ew^Si and ew^Ca contribution to the colour variability is performed sequentially, taking into account their correlation. Let the rest-frame colour ${c^{\prime }}_i(e{w}^{\mathit{Si}},e{w}^{\mathit{Ca}})$ be a function of the two equivalent widths with ${\sigma }_i=\partial {c^{\prime }}_i/\partial e{w}^{\mathit{Si}}$ and ${\kappa }_i=\partial {c^{\prime }}_i/\partial e{w}^{\mathit{Ca}}$. The small fluctuations (from SN to SN) are defined as δ ew = ew − ⟨ew⟩ (where ⟨ ⋅ ⟩ denotes the averaging over all SNe). The associated change in ‘colour difference’ $\delta c_i^{\prime }=e_i$ is given by the partial derivatives:

$$\begin{array}{c}\hfill e_i=\delta c_i^{\prime }={\sigma }_i\delta e{w}^{\mathit{Si}}+{\kappa }_i\delta e{w}^{\mathit{Ca}}.\end{array}$$(6)

As we are only interested in the derivatives, the subtraction of the averages is unnecessary. The coefficients σ_i and κ_i are not directly observed from the data as a consequence of the correlation between the two equivalent widths, which needs to be removed, and can be described by two linear relations:

$$\begin{array}{cc}& \delta e{w}^{\mathit{Si}}=a_{Si/Ca}\delta e{w}^{\mathit{Ca}}+\text{uncorrelated fluctuations,}\hfill \\ \hfill & \delta e{w}^{\mathit{Ca}}=a_{Ca/Si}\delta e{w}^{\mathit{Si}}+\text{uncorrelated fluctuations}.\hfill \end{array}$$(7)

The coefficients a_Si/Ca and a_Ca/Si are not inverse of each other as a consequence of measurement errors and uncorrelated physical fluctuation. They are found directly from the observed equivalent widths by two different linear fits to the observed distribution in the (ew^Si,ew^Ca) plane. Two χ² are minimised:

$$\begin{array}{cc}& {\chi }_{Ca/Si}^2={\displaystyle \sum _i{(e{w}_i^{\mathit{Ca}}-a_{Ca/Si}e{w}_i^{\mathit{Si}}+b_{Ca/Si})}^2,}\hfill \\ \hfill & {\chi }_{Si/Ca}^2={\displaystyle \sum _i{(e{w}_i^{\mathit{Si}}-a_{Si/Ca}e{w}_i^{\mathit{Ca}}+b_{Si/Ca})}^2.}\hfill \end{array}$$(8)

The dependence of ew^Si on ew^Ca gives a_Si/Ca = 0.0379 ± 0.0771, and the dependence of ew^Ca on ew^Si gives a_Ca/Si = 0.513 ± 0.286. The observed values are not weighted by the measurement errors, which do not include the intrinsic variability of the equivalent widths. The correlation a_Si/Ca of ew^Si as a function of ew^Ca is compatible with zero, and ignoring it would not have a significant impact on the corrections. We nevertheless take it into account to ensure that when the procedure described below is applied, there is no residual dependence whatsoever of colours on ew^Si or ew^Ca. We stress that we do not claim that either a_Si/Ca or a_Ca/Si represents the physical correlation between the equivalent widths, which is described by a single number $a_{Ca/Si}^{\text{phys}}=1/a_{Si/Ca}^{\text{phys}}$. Extracting $a_{Ca/Si}^{\text{phys}}$ would require a full understanding of the convolution of measurement errors and variability, which we do not actually need.

The large difference between the two coefficients reflects the large difference in the equivalent widths of the two lines. The observed dependencies of colours on the variations of Ca and Si equivalent widths are the total derivatives of the rest-frame colours with respect to ew^Si and ew^Ca, $K_i^{\mathit{Si}}$ and $K_i^{\mathit{Ca}}$, which can be expressed in terms of the partial derivatives ${\sigma }_i=\partial {c^{\prime }}_i/\partial e{w}^{\mathit{Si}}$ and ${\kappa }_i=\partial {c^{\prime }}_i/\partial e{w}^{\mathit{Ca}}$:

$$\begin{array}{cc}& \frac{d{c}_i^{\prime }}{de{w}^{\mathit{Si}}}=K_i^{\mathit{Si}}={\sigma }_i+{\kappa }_i{a}_{Ca/Si},\hfill \\ \hfill & \frac{d{c}_i^{\prime }}{de{w}^{\mathit{Ca}}}=K_i^{\mathit{Ca}}={\kappa }_i+{\sigma }_i{a}_{Si/Ca}.\hfill \end{array}$$(9)

The total derivatives K^Si and K^Ca are obtained from σ_i and κ_i in Eq. (9), once the coefficients a_Ca/Si and a_Si/Ca are known. Inversely, the coefficients σ_i and κ_i can then be derived in Eq. (10) from the observed dependencies of the colours on the equivalent widths, as seen in Figs. 1 and 2:

$$\begin{array}{cc}& {\sigma }_i=\frac{K_i^{\mathit{Si}}-a_{Ca/Si}{K}_i^{\mathit{Ca}}}{1-a_{Ca/Si}{a}_{Si/Ca}},\hfill \\ \hfill & {\kappa }_i=\frac{K_i^{\mathit{Ca}}-a_{Si/Ca}{K}_i^{\mathit{Si}}}{1-a_{Ca/Si}{a}_{Si/Ca}}.\hfill \end{array}$$(10)

Fig. 1. Correlations between the first five colours $c_0^{\prime }\equiv U-B$, $c_1^{\prime }\equiv U-V$, $c_2^{\prime }\equiv U-R$, $c_3^{\prime }\equiv U-I$, $c_4^{\prime }\equiv B-V$, $c_5^{\prime }\equiv B-R$ (mag) and ewSi (Å). The slopes of the straight line fits (red lines) give the value of KSi directly. The correlations of the colours can be derived from the first four.

Fig. 2. Correlations between the same colours as Fig. 1, corrected for the ewSi correlation (mag), as a function of ewCa (Å). The slopes of the straight line fits (red lines) give the value of K′Ca.

The denominator is zero when a_Ca/Si = 1/a_Si/Ca, but as explained previously, this would only occur if there were no fluctuations in the correlation between the two equivalent widths, and the issue addressed in this section would vanish.

We now use Eq. (10) to correct for the dependence of colours on the equivalent widths in the data. If the ew^Si-corrected colour is first applied, the observed dependence with ew^Ca will be ${K^{\prime }}_i^{\mathit{Ca}}={\kappa }_i(1-a_{Si/Ca}{a}_{Ca/Si})$. The values of K^Si and K′^Ca are directly obtained from the correlation between the colours and the equivalent widths observed in Figs. 1 and 2. The ew^Si and ew^Ca contributions to each colour $c_i^{\prime }$ are then subtracted to yield $c_i={c^{\prime }}_i-{\sigma }_{i}e{w}^{\mathit{Si}}-{\kappa }_{i}e{w}^{\mathit{Ca}}$. In all the following steps of the analysis, we use the colour differences e_i = c_i − ⟨c_i⟩, where ⟨ ⋅ ⟩ represents –as above– the mean over all the SNe Ia of our sample³. To preserve the algebraic relations between the ten colours, the ew^Si and ew^Ca corrections are performed on the first four colours, U − B, U − V, U − R, and U − I, while the others are obtained by combining them. The coefficients K^Si, K′^Ca, σ, and κ are summarised in Table 3 for the first four colours; the other colours can then be obtained by subtraction from the appropriate pair.

As shown in Sect. 9, we found that the correction of colours for spectral features was necessary. The parameters σ_i and κ_i could be considered as additional parameters of the model implemented, though they are not tuned but measured; we show in Sect. 9 that the Ca-Si colour correction is closely related to the SALT variable x₁.

4.2. Evaluation of extinction coefficients δ

4.2.1. From extinction formula to coefficients δ

This section describes the contribution of the reddening to the colour variations. Our reddening indicator is the extinction part X₄ (defined in Eqs. (3) and (4)) of e₄ = E(B − V)−⟨E(B − V)⟩ after the ew^Ca and ew^Si corrections have been applied; its mean ⟨X_i⟩ is zero, and it can have a positive or negative sign, but its variations are the same as the standard one (up to the intrinsic contributions of ew^Si, ew^Ca, and I_i, which we do want to exclude from the extinction correction). The SNfactory top-hat filters are defined in Table 1, and the extinction formula is (at first) assumed to be perfectly known, as given by FM99, once a value of R_V is chosen, which we assume until Sect. 7.

The reddening corrections $A_F^{\prime }$ and $R_{Fi}^{\prime }$ of Eq. (12) for the bandpass F are then evaluated for each value of R_V in a grid, and the analysis described in this section is repeated over the grid. The difficulty that we address is that the extinction indicator E_n(B − V) (with or without subtraction of the mean) is always tied to the choice of a wide bandpass, top-hat (our case), or Bessell. The reddening corrections at each wavelength depend on wavelength and R_V, but also on the choice of wide band filter, and on the spectral shapes of the selected stars through the wide-band integration. The rest-frame photon-count spectrum s_n(λ) for SN n is

$$\begin{array}{c}\hfill s_n(\lambda )=s_0(\lambda ){10}^{-0.4\varphi (\lambda )E_n(B-V)}.\end{array}$$(11)

With the colour difference definition E_n(B − V) = (B − V)_n − ⟨B − V⟩, the spectrum s₀(λ) at null colour difference is the mean SN Ia photon-count spectrum; the function ϕ(λ) is provided by Cardelli et al. (1989), O’Donnell (1994), Fitzpatrick & Massa (2005). This extinction indicator E_n(B − V), or ‘colour difference’, differs by a constant term from the standard E(B − V), which is referenced to the zero extinction value rather than the mean; the extinction indicator E_n(B − V) can be positive or negative. This constant shift (bandpass-dependent) affects colours, but not their variations, with which we deal exclusively. Equation (11) is a slight extension of the standard extinction formula of FM99 on two counts: the extinction formula relies on a specific set of blue stars, not SNe, and the bandpasses are defined by Bessell filters, rather than top-hat filters as described in Table 1. This issue is addressed in Sect. 4.2.2. Hereafter, we use the extinction component of E_n(B − V), that is, X₄(n) = X_B − V(n), as a substitute for the colour difference E_n in Eq. (11). It is shown in Appendix A that, in the linear approximation, the extinction correction $A_F^{\prime }$ in filter F is given for a SN with photon-count spectrum s(λ) by:

$$\begin{array}{cc}& A_F^{\prime }(n)=m_F^{\text{obs}}(n)-m_F^{\text{true}}(n)\hfill \\ \hfill & \simeq R_{F4}^{\prime }{X}_4(n)=\frac{{\int }_F{s}_0(\lambda )\varphi (\lambda )d\lambda }{{\int }_F{s}_0(\lambda )d\lambda }{X}_{B-V}(n).\hfill \end{array}$$(12)

The coefficient $R_{F4}^{\prime }$ is the exact derivative of the reddening correction $A_F^{\prime }$(n) with respect to the extinction colour component X₄ at X₄ = 0. We show in Table A.1 that the relative error resulting from this approximation remains smaller than 0.01 over the full range of values of X₄. The reddening correction coefficient $A_{F4}^{\prime }$ should actually be weighted by the spectrum of each SN according to Eq. (12), but as we only consider a linear reddening correction, the SN dependence would bring a negligible second-order correction. The term X_B − V(n) = X₄(n) in Eq. (12) is the extinction content of the colour difference e₄ = E_n(B − V) for SN n. As noted previously, the extinction colour X_i(n) can be positive or negative, but should be a better reddening indicator than E_n(B  −  V) as the intrinsic colour variation has been removed.

For each bandpass F = UBVRI, the extinction coefficient $R_{F4}^{\prime }=A_F^{\prime }/X_{B-V}(n)$ is obtained from Eq. (12). The associated colour differences are e_i = c_i − ⟨c_i⟩. The extinction coefficients $R_{Fj}^{\prime }$ (j ≠ 4) are then derived from the integrated extinction formula. The value of $R_{Fi}^{\prime }$ for another colour i can be derived from $R_{F4}^{\prime }$; for example $R_{F3}^{\prime }$ (with e₃ ≡ U − I) is given by:

$$\begin{array}{c}\hfill R_{F3}^{\prime }=\frac{R_{F4}^{\prime }}{R_{U4}^{\prime }-R_{I4}^{\prime }}.\end{array}$$(13)

The colour–colour coefficients ${\delta }_{ij}^{\prime }$ relate the (extinction) colour excesses of all colours; they are obtained from the two bandpasses defining each colour, from the analogue of Eq. (13); for example, for e₂ ≡ U − R and e₆ ≡ B − I, one has X₂ = δ₂₄X₄ = δ₂₆X₆ with

$$\begin{array}{c}\hfill {\delta }_{26}=R_{U6}^{\prime }-R_{R6}^{\prime }=\frac{R_{U4}^{\prime }-R_{R4}^{\prime }}{R_{B4}^{\prime }-R_{I4}^{\prime }}.\end{array}$$(14)

The extinction correction for all the other colours is similar. The colour indicators used in this study is actually the extinction component X_B − V of B − V computed in Sect. 5.2 rather than the measured value e_B − V. For instance, the correction to e₅ ≡ B − R using X₆ ≡ B − I is

$$\begin{array}{c}\hfill X_5={\delta }_{56}(R_V)X_6=\frac{R_{B4}^{\prime }-R_{R4}^{\prime }}{R_{B4}^{\prime }-R_{I4}^{\prime }}{X}_6={\delta }_{54}(R_V)X_4.\end{array}$$(15)

Using for instance the extinction part X₄ of colour 4 ≡ B − V or the extinction part X₆ of colour 6 ≡ B − I:

$$\begin{array}{cc}& m_F{(n)}^{\text{corr}}=m_F{(n)}^{\text{obs}}-A_F^{\prime }(n)(R_V)\hfill \end{array}$$(16)

$$\begin{array}{cc}& A_F^{\prime }(n,R_V)={\delta }_{F4}^{\prime }(R_V)X_4(n)={\delta }_{F6}^{\prime }(R_V)X_6(n).\hfill \end{array}$$(17)

(We recall that $m_F^{\text{obs}}(n)$ has been brought to a common reference redshift).

4.2.2. Rescaling of extinction coefficients: from δ′ to δ

The extinction formulae, such as those of Cardelli et al. (1989) or FM99, have not been tuned for SN spectra, nor to our set of top-hat bandpasses. In addition, the extinction law was established for a sample of mostly HD stars, where the extinction of each spectral bin is weighted by a spectrum that differs from the SNe Ia spectrum. It should not be assumed that the extinction formulae provided by Cardelli et al. (1989), O’Donnell (1994), Fitzpatrick & Massa (2005) –with different filters (Bessell) and different stars (HD)– apply directly to SNe Ia, given the presence of the wide band term E_n(B − V) or X₄(n) in Eq. (11). The minimal correction proposed here is a rescaling of the formula so as to ensure consistency. If we apply Eq. (12) to e₄ = E_n(B − V), we find $e_4=(R_{B4}^{\prime }-R_{V4}^{\prime })e_4$ (or $X_4=[R_{B4}^{\prime }-R_{V4}^{\prime }]X_4$). We must expect

$$\begin{array}{c}\hfill R_{B4}^{\prime }-R_{V4}^{\prime }=1.\end{array}$$(18)

If the extinction in the V bandpass is used –which is a frequent occurrence– a different consistency condition should be implemented: $R_{V4}^{\prime }=1$ to ensure that A_V is the actual extinction in this bandpass. Equation (18) is sometimes written (e.g. Mandel et al. 2017; Jha et al. 2007) as R_B − R_V = 1. However, we find (for R_V = 2.25) with the previously defined top-hat bandpasses ${(R_{B4}^{\prime }-R_{V4}^{\prime })}_{FM99}=0.864$, while with O’Donnell (1994) extinction parameters, $(R_{B4}^{\prime }-R_{V4}^{\prime })=0.814$. Most of this difference is due to the Bessell bandpasses assumed in the extinction formula, while we are using top-hat filters in the present study. Indeed, when using Bessell filter weighting (Bessell & Murphy 2012), $R_{B4}^{\prime }-R_{V4}^{\prime }=1.0329$, which is closer to but still incompatible with unity. We show in Eq. (A.5) of the Appendix that the expression of the coefficient $R_{F4}^{\prime }$ in Eq. (12) is the exact derivative of the reddening correction with respect to the extinction colour at X₄ = X_B − V = 0, and in Table A.1 that the linear approximation is accurate to within 5 10⁻³ mag over the full range of extinctions. Another contribution to this mismatch may be the different spectra of stars used in FM99 and in the present work (SNe Ia). We assume that an acceptable correction ensuring the ‘consistency condition’ in Eq. (18) is an overall rescaling of the extinction formula at all wavelengths, so as to force R_B4 − R_V4 = 1. All the coefficients $R_{Fi}^{\prime }$ are rescaled:

$$\begin{array}{c}\hfill R_{\mathit{Fi}}=R_{\mathit{Fi}}^{\prime }/(R_{B4}^{\prime }-R_{V4}^{\prime }).\end{array}$$(19)

The extinction coefficient of each filter is then increased by the factor $1/(R_{B4}^{\prime }-R_{V4}^{\prime })=1/0.864$. No such global rescaling was applied in C11 where an extinction parameter A_V was fit independently for each SN. As mentioned above, Eq. (11) is only an approximation. Even if one accepts that different molecular and grain compositions can be summarised by a single coefficient R_V (related to their cross section for light), its value could differ in the halo and in the disc (if there is one), and the observations of Schlafly et al. (2016) in our own galaxy suggest this possibility. The δ_ij coefficients of Eq. (4) verify that δ_ij = 1/δ_ji. The five values of the extinction coefficients δ_F4(R_V), which provide the extinction in filters for F = UBVRI from X₄ = X_B − V, are given in Table 4. The extinction coefficients δ_F(R_V) are subject to a first correction (Dex_F), which cancels the correlations of the colour residuals with the extinction components as discussed in Sect. 5.4; they are then rescaled for algebraic consistency as described in Eq. (19). Finally, for arbitrary values of R_V, the suppression of the correlation of extinction-corrected magnitudes with the extinction colour in the V bandpass requires an extra overall rescaling s, which is a function of the extinction parameter R_V, as explained in Sect. 7. Such a scaling destroys the consistency condition in Eq. (18), but is introduced in practice by almost all authors (Boone et al. 2021; Chotard et al. 2011; Amanullah et al. 2015; Saunders et al. 2018; Thorp et al. 2021) who substitute $m_F^{\text{corr}}(\lambda )=m_F^{\text{obs}}(\lambda )-A_V\varphi (\lambda ,R_V)$ where A_V is arbitrary and ϕ(λ, R_V) stands for one of the reddening formula. The analysis of the extinction is performed over a grid of nine values of R_V: 1.95, 2.05, 2.15, 2.20, 2.25, 2.35, 2.50, 2.60, and 3.10. We show below that for our value of R_V, the scaling factor s is indeed compatible with unity, which confirms the consistency of the model, as well as the value of R_V obtained.

5. Extraction of intrinsic and extinction colour components

5.1. Intrinsic colour coefficients γ

As mentioned in Eq. (3), the colour difference e_i(n) of SN n (after correction for ew^Si and ew^Ca) is the sum of two components e_i(n) = I_i(n)+X_i(n)+ϵ_i(n), and the intrinsic component I_i of colour i is assumed to belong to a one-dimensional space of intrinsic colour fluctuations, as suggested by Léget et al. (2020) and Boone et al. (2021). These are then interrelated by a set of coefficients γ_ij common to all supernovæ n as introduced in Eq. (4). As shown below, small corrections to the extinction formula are required. An initial value for the coefficients γ_ij can be derived directly from the requirement that intrinsic and extrinsic colours be uncorrelated, though this value can differ significantly from the result of the χ² minimisation discussed later. For any pair of leading (i) and auxiliary (j) colour, two equations link X_i, X_j and I_i, I_j:

$$\begin{array}{cc}\hfill e_i& =I_i+X_i+{ϵ}_i={\gamma }_{\mathit{ij}}{I}_j+{\delta }_{\mathit{ij}}{X}_j+{ϵ}_{\mathit{ij}},\hfill \end{array}$$(20)

$$\begin{array}{cc}\hfill e_j& =I_j+X_j+{ϵ}_j.\hfill \end{array}$$(21)

The initial values of γ_i0 can be obtained from the previous equations, assuming the absence of correlations between the I_i and the X_i, and neglecting the contribution of the residuals ϵ_i and ϵ_ij,

$$\begin{array}{cc}& ⟨e_1^2⟩={\gamma }_{10}^2⟨I_0^2⟩+{\delta }_{10}^2⟨X_0^2⟩,\hfill \end{array}$$(22)

$$\begin{array}{cc}& ⟨e_0^2⟩=⟨I_0^2⟩+⟨X_0^2⟩,\hfill \end{array}$$(23)

$$\begin{array}{cc}& ⟨e_0{e}_1⟩={\gamma }_{10}⟨I_0^2⟩+{\delta }_{10}⟨X_0^2⟩.\hfill \end{array}$$(24)

Here, δ_ij is the ratio of the extinction components X_i/X_j, as defined in Table 2 –and is not the Kronecker δ–, and the index n of the SN is dropped.

Combining Eqs. (23) (×δ₁₀) and (24), one gets $\langle e_0{e}_1\rangle -{\delta }_{10}\langle e_0^2\rangle =({\gamma }_{10}-{\delta }_{10})\langle I_0^2\rangle$. Similarly, ${\langle e_1\rangle }^2-{\delta }_{10}\langle e_0{e}_1\rangle ={\gamma }_{10}({\gamma }_{10}-{\delta }_{10})\langle I_0^2\rangle$; dividing the two previous relations one obtains an initial value of γ₁₀ (with large errors). More generally, an initial value for any pair of colours would be:

$$\begin{array}{c}\hfill {\gamma }_{\mathit{ij}}=\frac{⟨e_i^2⟩-{\delta }_{\mathit{ij}}⟨e_i{e}_j⟩}{⟨e_i{e}_j⟩-{\delta }_{\mathit{ij}}⟨e_j^2⟩}.\end{array}$$(25)

The initial values of the three coefficients γ₁₀, γ₂₀, and γ₃₀ are respectively 0.671, 0.759, and 0.401. The statistical errors on these initial values are relatively large, and an improved determination is obtained from the iterative fits described in the following section. (We recall that all colour differences have zero average, and that intrinsic and extinction components are assumed to be uncorrelated.)

Given the relation δ_ij = 1/δ_ji, we find that (as expected) γ_ij = 1/γ_ji. The coefficients γ_ij can be expressed according to Eq. (4) in terms of three of them, which have been selected as γ₁₀ = I_{(U − V)}/I_{(U − B)}, γ₂₀ = I_{(U − R)}/I_{(U − I)}, and γ₃₀ = I_{(U − I)}/I_{(U − B)}. The remaining 42 coefficients can be obtained from the following relations, which result from the algebraic relations between colours:

$$\begin{array}{cccccc}\hfill {\gamma }_{40}& =-1+{\gamma }_{10},\hfill & \hfill {\gamma }_{50}& =-1+{\gamma }_{20},\hfill & {\gamma }_{60}& =-1+{\gamma }_{30},\\ \hfill {\gamma }_{70}& =-{\gamma }_{10}+{\gamma }_{20},\hfill & \hfill {\gamma }_{80}& =-{\gamma }_{10}+{\gamma }_{30},\hfill & {\gamma }_{90}& =-{\gamma }_{20}+{\gamma }_{30},\\ & {\gamma }_{\mathit{ik}}\hfill & \hfill ={\gamma }_{i0}/{\gamma }_{k0}.\end{array}$$(26)

5.2. χ² for intrinsic and extinction colour components

The model described by Eq. (20) should reproduce the rest-frame colours. We consider the colour difference for SN n: e_i(n) = c_i(n)−⟨c_i⟩. The two components (I_i(n),X_i(n)) can be found by requiring that all colours e_j(n) are properly described by Eqs. (20) and (21), and the coordinates I_j(n) = γ_jiI_i(n), X_j(n) = δ_jiX_i(n). Colour components X₂(n) and I₂(n) are, for example, obtained by minimising –for each SN n– the ${\chi }_2^2(n)$ in Eq. (29) as a function of the ‘coordinates’ I₂(n) and X₂(n). The reddening correction A_V(n) in bandpass VA_V(n) = R_V4X₄(n), which varies from SN to SN, plays no role in the minimisation: the extinction correction is A_F4(n) = R_F4X₄(n) for filter F and SN n (=A_Fi(n) = R_FiX_i(n) using another colour):

$$\begin{array}{cc}& {ϵ}_2=X_2+I_2-e_2,\hfill \end{array}$$(27)

$$\begin{array}{cc}& {ϵ}_{k2}=({\gamma }_{k2}{I}_2+{\delta }_{k2}{X}_2-e_k),\hfill \end{array}$$(28)

$$\begin{array}{cc}& {\chi }_2^2={ϵ}_2^2+{\displaystyle \sum _{k\in Aux}{ϵ}_{k2}^2,}\hfill \end{array}$$(29)

as a function of I₂ and X₂, under the constraint

$$\begin{array}{c}\hfill {\displaystyle \sum _{\text{all SNe}}{I}_2(n)X_2(n)=0.}\end{array}$$(30)

Equation (30) imposes the decorrelation of the extra-intrinsic colour from extinction. The index n is omitted in Eq. (27). For each colour i and each SN, the residual ϵ_i(n) defined by Eq. (27) should not be confused with the magnitude shift at maximum light defined in Eq. (2). We are aware that Eq. (30) could be criticised under certain circumstances, such as in the presence of circumstellar matter with a composition differing from the average host galaxy. Such matter, as pointed out by Borkowski et al. (2009), Ferretti et al. (2017), might invalidate this assumption, but if the circumstellar dust is created by a non-degenerate companion, as suggested in one case by Nagao et al. (2017), it is not necessarily correlated to the SN Ia properties. While not rejecting the possibility, we wanted to explore the accuracy that could be reached within the assumption of decorrelation.

Each pair (I_i(n),X_i(n)) of components allows the prediction of all the other colours of the SN. When defining χ²(n) in Eq. (29), each colour i enters with equal weight, which implicitly assumes that the errors on all colours are similar. The solution of Eq. (29) in terms of (X₂(n),Y₂(n)) is found by requiring $\partial {\chi }_2^2/\partial X_2=\partial {\chi }_2^2/\partial I_2=0$:

$$\begin{array}{c}\hfill \left(\begin{array}{c}{Y}_0\\ Y_1\end{array}\right)=\left(\begin{array}{cc}{M}_{00}& M_{01}\\ M_{10}& M_{11}\end{array}\right)\left(\begin{array}{c}{I}_2\\ X_2\end{array}\right)\end{array}$$(31)

with

$$\begin{array}{cc}& M_{00}=1+{\displaystyle \sum _{k\in Aux}{\gamma }_{k2}^2,M_{01}=1+\sum _{k\in Aux}{\gamma }_{k2}{\delta }_{k2}-\lambda ,}\hfill \\ \hfill & M_{10}=M_{01},M_{11}=1+{\displaystyle \sum _{k\in Aux}{\delta }_{k2}^2,}\hfill \\ & Y_0=e_2(k)+{\displaystyle \sum _{k\in Aux}{\gamma }_{k2}{e}_k,Y_1=e_2(k)+\sum _{k\in Aux}{\delta }_{k2}{e}_k.}\hfill \end{array}$$(32)

The colours excesses e₁(n) and e₂(n) are measured, meaning that Y₀ and Y₁ are known. The intrinsic and extinction components I₂(n), X₂(n) are obtained by inverting Eq. (31). The Lagrange multiplier λ arises from the assumption of uncorrelated I_i, X_i, and is found from Eq. (30); it depends on the colour, but is the same for all SNe. We introduce the set Aux of auxiliary colours used to extract the intrinsic and extinction components of each colour. In the optimisation of the intrinsic couplings γ_ij below, the full set of ten colours cannot be used in view of the appearance of algebraic constraints, as explained in the following section. Table 5 presents the actual list of auxiliary colours used for each colour. There is an arbitrariness in the choice of the participating colours in Aux; we selected the one that yielded the smallest residuals among five trials.

For each SN n, this minimisation yields ten intrinsic colours I_i(n), and ten extinction colours X_i(n), once the intrinsic colour coefficients (γ₀₁, γ₀₂, γ₀₃) are known. Equation (29) does not take into account the correlation between colours, but as the solution gave acceptable residuals for all colours, we kept this simplified expression for χ². We verified that including the obvious algebraic correlations between the different colours does not change the solution, but it does add significant complications. As the mean of e_i is zero (by definition), the mean values of I_i and X_i obtained by linear combinations of e_i will also be zero.

5.3. Determination of intrinsic components and intrinsic couplings

The colour differences which we now use were subtracted from the contributions of Ca and Si, but leave room for an extra-intrinsic component I_i of each colour, with all of them dependent on three parameters (γ₁₀, γ₂₀, γ₃₀) as required by the assumption of a single extra-intrinsic contribution. Up to this stage, the γ_ij are obtained from their initial values in Eq. (25). The solution of Eq. (31) is only acceptable if the ratio of the solutions of Eq. (31) I_i(n) to I_j(n) for two colours i and j is equal to γ_ij for all SNe and all pairs of colours. Similarly, the ratio of X_i(n) to X_j(n) should be equal to δ_ij.

This is achieved by a sequence of iterations. At each iteration, and for each colour, the input coefficients ${\gamma }_{\mathit{ij}}^{\text{in}}$ are the parameters introduced in Eqs. (4) and (28), while the linear coefficients of the straight line fits in the plane (I_i, I_j) as in Fig. 3 are the output value ${\gamma }_{\mathit{ij}}^{\text{out}}$ of the intrinsic couplings. The output values of (γ₁₀, γ₂₀, γ₃₀) are reinjected into Eq. (29) until Fig. 3 is satisfactory. The scatter is seen to range from 0.002 mag to 0.01 mag depending on the colours, and the measured values of the ratios ${\gamma }_{\mathit{ij}}^{\text{out}}=I_i/I_j$ are within 0.01 of the corresponding set of ${\gamma }_{\mathit{ij}}^{\text{in}}$, as obtained from (γ₁₀, γ₂₀, γ₃₀) using Eq. (26). The δ_ij coefficients are left unchanged in this sequence, as their values are fixed by the extinction formula (at this stage).

Fig. 3. Intrinsic ${\gamma }_{0i}^{\text{out}}=I_0/I_i$ and extinction ${\delta }_{0i}^{\text{out}}=X_0/X_i$ ratios as defined in Eq. (31) (i = 2, 6, 9) compared to the input ${\gamma }_{0i}^{\mathrm{in}}$ and calculated δ0i for RV = 2.25. The compatibility is at the level of 0.001 for the extinction, and the relative accuracy of 0.05 for the (smaller) intrinsic colour (worst case). The red lines are linear fits to the data, and γout, δout are the values of the linear coefficients.

We then turn to Fig. 4, which compares the full set of input values ${\gamma }_{\mathit{ij}}^{\text{in}}$ in Eq. (31) in the second-to-last iteration to the output values, that is, of I_i/I_j from the linear fits in Fig. 3. The final values of the three coefficients (${\gamma }_{10}^{\text{out}}$, ${\gamma }_{20}^{\text{out}}$, ${\gamma }_{30}^{\text{out}}$) are found by minimising the scatter (RMS) of the output coefficients $\sigma ({\gamma }_{\mathit{ij}}^{\text{out}})$ with respect to the input coefficients (horizontal scale). This last step brings in the constraint of consistency with Eq. (26) for the whole set of γ_ij, and modifies the result of the previous sequence of iterations; however, as seen in Fig. 3, the outcome is still satisfactory with respect to the convergence of the solution for each γ_ij

Fig. 4. Relation between input and output (⟨Ii/Ij⟩) values for the coefficients γij. The overall agreement for all γij (excluding those involving colours 5 ≡ (B − R) and 7 ≡ (V − R)) is within 0.01.

The relation between ${\gamma }_{\mathit{ij}}^{\text{out}}$ and ${\gamma }_{\mathit{ij}}^{\text{in}}$ is shown in Fig. 4 at the last iteration. The linear coefficient is expected to be close to unity if convergence has been achieved, while we find 1.00069, with an RMS of 0.0104. This dispersion arises from a combination of measurement errors, Si-Ca correction errors, and modelling errors. Two precautions were implemented to extract the three intrinsic couplings from Eq. (31) and Fig. 4: we eliminated the auxiliary colours 4 ≡ B − V and 7 ≡ V − R that have very small intrinsic components (cf. Table 9) from the 10 × 9/2 = 45 γ_ij couplings, 28 remain after the elimination of colours 5 and 7. The 28 output values of the parameters ${\gamma }_{\mathit{ij}}^{\text{out}}=I_i/I_j$ at the last iteration are compared in Fig. 4 to the input values derived from (γ₁₀, γ₂₀, γ₃₀) using Eq. (26). The equality of input and output values was not forced into the optimisation process. The ratios δ_ij = X_i/X_j are also seen to have a small scatter of 0.002 mag (in Fig. 3), with averaged output values close to the input coefficients, supporting the model Eqs. (27)–(29).

Nevertheless, the RMS in Fig. 4 fails as an indicator when the same set of colours (auxiliary and main) is used to find the intrinsic content, as the algorithm defined by Eq. (31) yields $\sigma ({\gamma }_{\mathit{ij}}^{\text{out}})=0$. For each colour, we eliminated at least one extra auxiliary colour to avoid such duplicated sets. The choice of the auxiliary colours has an impact on the values found for the three intrinsic couplings. We selected the one that gave the smallest mean residual scatter (RMS) in Fig. 4. The scatter of ${\gamma }_{\mathit{ij}}^{\text{out}}$ in Fig. 4 is 0.0103, and the ratio ${\gamma }_{\mathit{ij}}^{\text{out}}/{\gamma }_{\mathit{ij}}^{\text{in}}$ is within 0.001 of unity. After the sequence of iterations described in Sect. 5.4, the optimal values in Fig. 4 for R_V = 2.20 are given in Table 6.

The errors are derived at this stage from the spread of the results in four different choices of auxiliary colours and do not include the statistical error. The comparison of our intrinsic colours to previous results is postponed to Sect. 9, as the contribution of the equivalent widths ew^Ca and ew^Si must be reintroduced.

5.4. Modifications to the extinction formula

For all values of R_V, the colour residuals ϵ_i described in the previous section (after Eq. (29)) show a strong correlation with the extinction colour component when the FM99 extinction formula is used. This (unexpected) correlation is apparent in Fig. 5, and there is no value of R_V that suppresses this effect. As our subsequent determination of R_V relies on the absence of any correlation between the extinction-corrected magnitudes and the extinction, it is crucial to suppress the correlation observed in Fig. 5. This suppression depends on the coefficients γ_ij and δ_ij in Eq. (32) but in no way on the scale of the reddening correction (nor on the rescaling applied in Eq. (19)).

Fig. 5. Correlation of the colour residuals (mag) for RV = 2.25 with the extinction colour for colour differences e0, e2, e5, e6, e8, e9 (no extinction formula correction). The strong correlations observed would forbid a safe evaluation of RV.

We can cancel the effects observed in Fig. 5 by introducing corrections Dex_F to the extinction coefficients R_F4 of the five F = UBVRI bandpasses in Eq. (12). The minimised quantity is

$$\begin{array}{c}\hfill {\chi }_{\mathrm{Dex}}^2={\displaystyle \sum _{i=0}^9{(s{l}_i)}^2\text{with}s{l}_i=\partial {ϵ}_i/\partial X_i.}\end{array}$$(33)

The extinction corrections Dex_F in the five filters are applied to the values of $R_{F4}^{\prime }$ from Eq. (12): $R_{F4}^{\prime }\to R_{F4}^{\prime }+De{x}_F$ and fitted for each value of R_V so as to minimise ${\chi }_{Dex}^2$, and cancel the coefficients sl_i in Fig. 5. The result is given in Table 7. As there is an overall scale degeneracy of the extinctions in the colour analysis, we arbitrarily set the U bandpass correction to be 0.0675. At each value R_V of the grid, the corrections are evaluated. As these corrections (partially) compensate for the distortions introduced by the difference between the ‘optimal R_V’ and the value used in the grid, they tend to blow up when R_V is close to the limits of the range considered in this work, namely as R_V → 1.95 or R_V → 3.10. The offsets are applied in the second line of Table 4. The impact of these corrections on the correlation between the residuals and the extinction is shown in Fig. 6, and for all colours in Table 8. The large value of sl₉ observed (e₉ ≡ R − I) in Fig. 5 reflects the significant corrections to the extinction formula in bandpasses R and I, with opposite signs. The errors given in these figures are only indicative, as somewhat arbitrary ‘floor’ errors of 0.006 mag are imposed on the determination of the extinction colour.

Fig. 6. Correlation of residuals of colour differences (mag) (e0, e2, e5, e6, e8, e9) with extinction after correcting the extinction formula with RV = 2.25. There are only four independent colours, meaning that the four extinction corrections allow us to cancel the correlations of all ten colours.

The corrections to the extinction formula are extremely unlikely to be caused by the limitations of our model: The linear approximation involved is numerically correct to better than 0.005 mag in the whole range of e₄ = E(B − V) considered. There could be a second intrinsic colour involved; that is, rather than one-dimensional (beyond Ca and Si), as in Eq. (4), the intrinsic colour space could be two-dimensional, with a second set of γ_ij coefficients, but intrinsic colours are uncorrelated (or are at best weakly correlated) to extinction, and this new intrinsic colour would have to be larger than the one already introduced to make up for the 0.10 mag discrepancies observed in Fig. 5. There is no room in the residuals for such a large extra-intrinsic component. On the other hand, the need for an adaptation of the extinction formula to SNe and Top-hat filters is shown in Sect. 4.2.2, where the presence of a rescaling factor is first introduced.

The numerical value of the corrections are dependent on the specific bandpasses involved, on the spectral template, and also possibly on the dust properties of our galaxy, which may differ from the averaged dust extinction. Figure 6 shows that, after correction, the remaining correlations of the residuals ϵ_i with extinction are negligible. This can be achieved regardless of the value of R_V. The largest of these effects amounts to a residual linear dependence of 0.006 for B − V. The coefficients of all other colours are smaller than 0.0036. The statistical errors on the offsets away from the extinction formula for our sample of SNe will be estimated from a simulation in Sect. 8.4.

It seems unnatural to constrain the γ_ij coefficients with ${\chi }_{Dex}^2$, although they are (weakly) sensitive to the value of the extinction corrections. Our choice is instead to optimise the dispersion σ(γ) displayed in Fig. 4, as stated in Sect. 5.3. At each value of the parameters (Dex_F or γ_i0), the I_i and X_i are obtained from Eq. (31), and the γ_i0 are retuned in order to minimise σ(γ) in Fig. 4. This sequence is repeated four or five times to reach the final values of the seven unknown parameters (four Dex_F and three γ_i0). Within each cycle (extinction corrections Dex_F or γ_i0 fixed), convergence is reached when ${\chi }_{Dex}^2$ or σ(γ) increase whenever any one of the seven parameters is changed by 10⁻³. As the implemented model is still imperfect, the minima of ${\chi }_{Dex}^2$ and σ(γ) are not perfectly compatible; the differences between the two minima are nevertheless approximately ten times smaller than the quoted errors.

6. Quality of the colour reconstruction

The quality of the reconstruction of the colour differences e_i from its intrinsic I_i and extinction X_i components is shown in Table 9 and in Fig. 7, where the distribution of the difference ϵ_i = c_i − (I_i + X_i) is given for six selected colours. The average colour residual over the ten colours is 0.0176. As the colour measurement errors do not exceed 0.006, this figure (together with Table 9) suggests that the dominant contributions arise from modelling error and from the subtraction of the ew^Si and ew^Ca contributions. Given the accuracy reached, the phase of the spectrum might also contribute.

Fig. 7. Residuals ϵi (mag) after the reconstruction of six colours (c0, c2, c5, c6, c8, c9). The full list of residuals is given in Table 8.

7. Determination of R_V

7.1. Reddening correction of magnitudes

In Sect. 5.4, we ensure that the colour residuals are independent of the extinction component for each selected value of R_V. However, after correction for reddening, the magnitudes and the colours themselves may still depend on the extinction –even though the dependences in data and model must be the same after the corrections of Sect. 5.4. For each choice of R_V, the magnitude m_F found from the data at a reference redshift of 0.05 is corrected for its dependence on ew^Si and ew^Ca –as carried out in Sect. 4.1 for colours– in order to obtain $m_F^{\mathit{CaSi}}=m_F-{\alpha }_{\mathit{Si}}e{w}^{\mathit{Si}}-{\beta }_{\mathrm{Ca}}e{w}^{\mathit{Ca}}$. The dependence of $m_F^{\mathit{CaSi}}$ on the intrinsic component I_i with the (measured) coefficient z_F1 is then taken into account, and we allow for a single overall scaling factor s(R_V), which mimics the arbitrary A_V used in most analyses (it should be unity in our framework) to define the extinction-corrected magnitude of each SN in the bandpass F using colour i:

$$\begin{array}{c}\hfill m_F^{\text{corr}}=m_F^{\mathit{CaSi}}-z_{F1}{I}_1-s(R_V)R_{F1}(R_V)X_1.\end{array}$$(34)

Once the extinctions have been scaled according to Eq. (19), as explained in Sect. 4.2.2, the magnitude $m_F^{\text{corr}}$ should not depend on X_i. The extinction scale s(R_V) in Eq. (34) is chosen so as to cancel the dependence of the corrected data magnitude in the V (top-hat) bandpass, and we expect s(R_V) = 1 for the physical value of R_V. It is seen in Fig. 8 that this is true when R_V ≈ 2.25. The extinction in filter F is now A_F = s(R_V)×R_F3 × X₃. However, if we turn to the other filters, we see in Fig. 9 that even after the extra scaling factor s(λ) has been introduced in Eq. (34) to ensure $d{\text{mag}}_V^{\text{corr}}/d{X}_3=0$, the derivatives of the other bandpasses are finite when R_V lies away from the 2.15 − 2.3 range. The dependence of $d{m}_F^{\text{corr}}/d{X}_3$ as a function of R_V is shown for the bandpasses UBRI in Fig. 9. We find the extinction scale s to be s(R_V = 2.20) = 1.0269 and s(R_V = 2.25) = 1.0014 (close to our ‘effective’ value of R_V). The errors in Figs. 8–11 are strongly correlated over the range of wavelengths, and for different values of R_V as the same events are used. The error is evaluated from the dispersion of the scale factor in the ten samples of the simulation. From sample to sample, the scale parameter moves up and down in Fig. 8, but the smooth behaviour is preserved as a function of R_V as the same events are involved for a given sample.

Fig. 8. Extinction scale factor s as a function of RV for dmV/dX3 = 0. For each value of RV, there is a value of the extinction scale factor s that cancels the dependence of the V magnitude on extinction.

Fig. 9. $d{m}_F^{\text{corr}}/d{X}_3$ for the UBVRI bandpasses as a function of RV. Although the dependence of the V bandpasses as a function of extinction has been cancelled for each value of RV, the extinction in other bandpasses is not adequately corrected.

Fig. 10. d(U−I)corr/dX3 as a function of RV. The difference between U and I magnitudes (colour 3) is the most sensitive to the choice of RV. The error bar is the dispersion found in the simulation.

Fig. 11. UBVRI magnitudes corrected for I3, X3 with RV = 3.100, s = 0.73344. The derivative $d{\text{mag}}_F^{\text{corr}}/d{X}_3$ is compatible with zero in each bandpass, but the variation from U to I is quite significant, taking correlations into account. The red line is weighted by errors and the blue straight line is unweighted.

Each bandpass provides a value of the parameter R_V that cancels the derivative $d{m}_F^{\text{corr}}/d{X}_3$. For UBRI, these values are respectively 2.283, 2.183, 2.137, and 2.240. We average the U and I measurements –which have a stronger R_V dependence– and obtain (with the error from the simulation):

$$\begin{array}{c}\hfill R_V=2.262\pm 0.126.\end{array}$$(35)

The scale factor s for this value of R_V is s = 1.001 ± 0.045, which is additional independent confirmation of the result for R_V. We used magnitudes rather than colours to derive the previous result, but it can be reframed as a colour measurement: Fig. 9 shows that the largest mismatch in the value of $d{m}_F^{\text{corr}}/d{X}_3$ is obtained by comparing the U and I bandpasses. We take advantage of this observation in Fig. 10 by imposing the condition that the derivative of the U − I colour with respect to the extinction colour X₃ be zero after application of the extinction correction. The value of R_V obtained in Fig. 10 is now R_V = 2.265, which is almost the same as in Eq. (35), but the error is slightly smaller, as expected; the contributions from flux calibration error (0.03 mag), the redshift error, and the ‘grey’ fluctuation are suppressed. With this determination, the corresponding scale factor for the ‘extinction scale’ s = 1.000. If we turn to the usual ‘operational’ definition of R_V, as A_V/E(B − V), it is seen in Table 4 that R_V = A_V/X₄ = R_V4 = 2.7. However, this numerical value is specific to SNe and our top-hat filters.

The corrected magnitudes in the UBVRI bandpasses are shown in Fig. 12 for the value R_V = 2.25. As expected, no residual dependence of the (corrected) magnitude on the extinction is observed in any filter. The RMS (which was not minimised), has a value of 0.13 mag in all filters, which is smaller than the scatter of 0.15 mag typically seen in the SALT2 analyses of our sample (Saunders et al. 2018), but larger than the result from refined analyses designed to minimise this fluctuation, as in Boone et al. (2021). The blue straight line in Fig. 12 is a weighted fit and the red line is unweighted; the two slopes are almost identical. The magnitude fluctuation is ‘grey’, that is, it is almost identical over all bandpasses to a remarkable accuracy. For comparison, we show the magnitude correlations obtained for R_V = 3.1 in Fig. 11. The scale factor must now be set to s = 0.857 to ensure that the derivative $d{m}_V^{\text{corr}}/d{X}_3=0$. The error on the derivative of d(m_U−m_I)/dX₃ in Fig. 10 is 0.04, as evaluated from the spread of the simulation samples, meaning that the significance of the rejection of R_V = 3.1 (and higher) is actually 3.5 standard deviations.

Fig. 12. UBVRI magnitudes corrected for I3, X3 with RV = 2.25, s = 1.0024. As for the previous figure, the red straight line is weighted by errors.

7.2. Deriving R_V from the grey fluctuation

We mention in the previous section that the fluctuation of the magnitudes around the common ‘grey’ fluctuation of each SN Ia i, $m_{\text{grey}}(i)=\left(m_U^{\text{corr}}(i)+m_B^{\text{corr}}(i)+m_V^{\text{corr}}(i)+m_R^{\text{corr}}(i)+m_I^{\text{corr}}(i)\right)/5$, was remarkably small. This observation can be turned around and used to derive R_V: in each bandpass F, we compute the RMS σ_F of the dispersion of the corrected magnitudes around m_grey. We show how this dispersion varies with R_V in Fig. 13: a minimum is reached for R_V = 2.277, which confirms the previous result. The dispersion in each UBVRI bandpass with respect to the average offset over all filters is remarkably small, that is 0.0114, 0.0197, 0.0145, 0.0148, and 0.0060, respectively, which is everywhere smaller than 0.02. This result quantifies the level at which the remaining magnitude fluctuation is grey, and also confirms the validity of the modelling implemented (up to this ‘unaccounted’ grey fluctuation).

Fig. 13. RMS of the offsets (with respect to mgrey) of the corrected UBVRI magnitudes as a function of RV.

8. Simulation

The previous results can be summarised by ten numbers: the five extinction corrections to the extinction formula in the five bandpasses used for our analysis, the three intrinsic colour ratios, the extinction scale s, and the extinction parameter R_V. These are the ‘effective’ values found within our algorithms, but not necessarily the ‘true’ values. The biases and the errors are evaluated below based on a simulation that mimics the observations. We use the values of the parameters tuned to describe the data to generate the events, and then proceed with the same analysis. The differences from the reconstructed values allow us to find the methodological biases, and the scatter among the different generated samples helps us find the statistical error.

8.1. Generation of colours with their error

We generate nine samples of 165 SNe Ia UBVRI magnitudes. The generated intrinsic $I_1^g$ and extinction $X_1^g$ components of colour 1 ≡ U − V follow the observed distribution in the data after smoothing by a sliding local averaging (to avoid counting the statistical fluctuations twice): in each histogram bin, the data distribution of I₁ and X₁ was replaced by the average of three adjacent bins. Whenever the bin statistics was below four events, the average of five bins was used. The other intrinsic colours $I_k^g$ obey $I_k^g={\gamma }_{k1}{I}_1^g$, and the extinction colours $X_k^g$ follow from $X_k^g={\delta }_{k1}{X}_1^g$. The coefficients δ_ij used in the generation are the ones found in the data analysis; that is, they include the correction to the FM99 formula described in Table 7. The simulation shows that the analysis does recover the input corrections, with a small bias. As the measurement errors are included in the observed data distribution, they should have been removed from the generation in the simulation. We instead introduced a single simulation scaling factor of the $I_1^g$ distribution, tuned so that the reconstructed intrinsic components should match the observations once the analysis is applied; it depends on the simulated sample, but its difference from unity never exceeds 0.01. This effect is expected to be smaller for the extinction colours X_i, which have a larger range, and the observed distribution of the data has been used.

8.2. Colour noise

The noise is derived from the residuals observed in the data, and it combines measurement and modelling errors. The residuals ϵ_i measured from the data in four colours (i.e. 2 ≡ U − R, 5 ≡ B − R, 7 ≡ V − R, 9 ≡ R − I) are used to derive a 4 × 4 covariance matrix ⟨ϵ_iϵ_j⟩. The square roots of the dominant eigenvalues are σ₁ = 0.026525, σ₂ = 0.017947; the other eigenvalues are negligible (∼10⁻⁴ and ∼3 10⁻⁵). The Gaussian noise of the two eigenvectors is then projected onto the four bandpasses to obtain a simulated ‘colour noise’ for the four colours. We find that this noise has to be multiplied by 1.004 in the simulation to reproduce the observed residuals. This effect may arise from the contribution of modelling errors to the data, while the model is exact in the simulation.

8.3. Generation of magnitudes

The R magnitude is used as the reference. A correction reproducing its (small) correlation with the intrinsic colour is added. The magnitudes in the UBVI bandpasses are obtained by adding the U − R, B − R, V − R, and I − R colours, as found from their intrinsic and extinction components, and including the ‘colour noise’. Finally, a Gaussian grey magnitude fluctuation with σ(grey) = 0.12 is added, though its value is irrelevant for all the results of the present study. The generated magnitudes of event n in bandpass R and U are therefore described as

$$\begin{array}{cc}& m_R(n)=z_{R1}{I}_1(n)+R_{R1}{X}_1(n)+\text{grey}(n),\hfill \end{array}$$(36)

$$\begin{array}{cc}& m_U(n)=m_R(n)+{\gamma }_{21}{I}_1(n)+{\delta }_{21}{X}_1(n)+{\sigma }_1{u}_1(n)n_{21}\hfill \\ \hfill & +{\sigma }_2{u}_2(n)n_{22},\hfill \end{array}$$(37)

where u₁(n) and u₂(n) are centred random Gaussian variables with a standard deviation of unity, n₂₁ and n₂₂ are the components of the corresponding eigenvectors on the U − R colour, and z_R1 is the observed correlation between the R magnitude and the intrinsic colour 1. The generated magnitudes are processed by the same algorithms as the data (after ew^Ca and ew^Si corrections), providing reconstructed values for I_i and X_i, reconstructed extinction corrections, reconstructed intrinsic couplings, and a reconstructed R_V.

8.4. Simulated and observed distributions

As we want to evaluate the errors and biases of the data analysis, with the help of the simulation, both the distributions and the errors must be similar in the observations and the simulation. The (I_i, X_i) obtained through the colour analysis (as in the data), and the generated $(I_i^g,X_i^g)$ values are compared for three colours in Fig. 14 for the simulated sample 1 and in Table 9. There is a very small scatter between the generated and reconstructed colour components, and their ratio is close to unity. This scatter varies from 0.002 to 0.006 for I_i, and from 0.001 to 0.006 for X_i. The systematic trend in the ratio of measured and generated colour is partially statistical and fluctuates from one sample to another in a range of 0.02, which justifies using the observed distributions as an input for the simulation. The scatter in Fig. 14 is surprisingly smaller than the RMS of the observed residuals (in average 0.0176 mag). This is a consequence of the averaging over different colours performed in Eq. (31). The simulated accuracy is tuned to reproduce the observed residuals by the introduction of two ad hoc factors in the noise and the rescaling of the distribution of I₁ (of colour U − V) used in the generation. The distribution of the residuals is shown in Fig. 15 and Table 9.

Fig. 14. Comparison of intrinsic and extinction components (mag) in generation and reconstruction for colours c0, c5, and c9. The small differences found justify using the (smoothed) observed distributions of I1 and X1 as input.

For the same three colours (c₀, c₅, c₉), the intrinsic and extinction colour distributions for the data and one of the simulated samples are shown in Fig. 16. The results for all colours are summarised in Table 9. The residuals in the data and the simulation are compared in Fig. 15 and Table 9. The intrinsic and extinction colours in the simulation and observation agree, and the trends of the errors are also well reproduced by the simulation.

Fig. 15. Residuals in observations and simulations for colours c0, c2, c5, c6, c8, and c9 (RV = 2.25). The RMS of the residuals of all colours are given in Table 9. The blue shade is the simulated distribution, the orange shade the data.

Fig. 16. Distributions of intrinsic and extinction components of c0, c6, and c9 in observation and simulation for RV = 2.25 (mag). The extinction is asymmetric, with an extended tail to high values. The intrinsic part is symmetric. The RMS for all colours is given in Table 9. The blue shade is the simulated distribution, the orange shade the data.

8.5. Errors on the extinction formula corrections

We rely on the nine simulated samples to control the bias in the measurements arising from the implemented algorithms, as well as the errors on the extracted parameters: extinction corrections Dex_F, the value of R_V, and the extinction scale s_X. The analysis was performed on each of the samples, on the same grid of values of R_V. The mean value of the reconstructed extinction corrections and their dispersion is given in Table 10.

When comparing to Table 4, it is seen that the corrections to the extinction are statistically significant, rising from 1 to 10 percent over the visible spectrum once the bias is taken into account. The algorithm used in deriving the corrections is able to extract their value, with some biases. Accepting the biases found in the simulation, our final values for the corrections to the extinction coefficients are given in the last column of Table 10 for R_V = 2.25. The error on the bias correction is σ(Dex_F)/3 as nine samples are generated.

8.6. Errors on the intrinsic coupling coefficients γ

For each simulation sample and each value of R_V in the grid, the three intrinsic coupling coefficients γ₁₀, γ₂₀, and γ₃₀, are evaluated as in the data analysis, that is by minimising the scatter in the plot comparing the input values to the results. The standard deviation of the three coefficients from one sample to the next is respectively 0.0075, 0.0061, and 0.0041. The values averaged over all samples allow us to estimate the bias arising from the analysis; it is found to be of the same order as the statistical error. The largest contribution to the error is by far the one discussed in Sect. 5.1, which arises from different choices of auxiliary colours, and is given in Table 6. The statistical error will be added into the final uncertainty.

8.7. Errors on the extinction parameter R_V

The measurement of R_V is also performed on the simulation. The value in the generation was R_V = 2.20, and the average value found is R_V = 2.284 ± 0.112. We correct the result obtained in the data analysis, of 2.265, to R_V = 2.265 − 0.084 = 2.181 ± 0.117 (the error is increased to account for the uncertainty on the bias correction). As the values found in the four filters BVRI are almost fully correlated, we use the error on R_V derived from the U − I colour. We add the impact of an estimated systematic error in the measurement of colours and other corrections (earth atmosphere, galactic extinction, host galaxy subtraction, silicon and calcium contribution). As we only use differences between mean and measured colours, many instrumental systematic errors cancel out. To evaluate the impact on R_V, all colours were modified by a shift proportional to their value and reaching 0.005 mag over their range of variation. The outcome is a change in R_V by 0.05. The final result for the parameter R_V (including the error on the bias) is then

$$\begin{array}{c}\hfill R_V=2.181\pm 0.117(\text{stat.})\pm 0.050(\text{syst.}).\end{array}$$(38)

9. No calcium or silicon line corrections

We tried to reproduce our analysis without performing the initial correction for the ew^Si and ew^Ca spectral lines. The intrinsic colours are expected to be much larger than in the previous analysis. They reach 0.10 in colours (U − B, U − V, U − R). The mean colour residual is then 0.0402 (instead of 0.0176 with the colour corrections from ew^Ca and ew^Si). The square roots of the dominant eigenvalues of the covariance matrix from the residuals are 0.0724 and 0.0296, meaning that it would be necessary to add two spectral coordinates for each SN in order to reduce these residuals, instead of tolerating these contributions as ‘noise’ as is done in the present work. The two eigenvectors are (0.5842, 0.5974, 0.4373, −0.3324) and (0.5972, −0.7405, 0.2908, 0.1013). In the absence of ew^Si and ew^Ca corrections, the model we used fails to describe the observations, as we are missing key spectral information. Instead, the approach we take is to add back the colour contributions of Si IIλ4131 and Ca II H&Kλ3945 –which are subtracted in Sect. 4.1– to obtain a ‘full’ intrinsic component of colour variability; this is shown in Fig. 17.

Fig. 17. Distributions of six ‘full’ intrinsic colours including Si II and Ca II H&K contributions.

The range of the full intrinsic colour component is given in Table 11 for all colours. The small size of the intrinsic part of e₄ ≡ B − V justifies using it as a simplified indicator of extinction for SNe Ia.

The standard choice of colour 4 ≡ B − V as an extinction indicator is fortunate, as the intrinsic content happens to be small. However, a large part of the calcium and silicon correction can be recovered, with information extracted from the light curve. The SALT2 stretch variable x₁, as described in Betoule et al. (2014), is strongly correlated to the combined calcium–silicon corrections DCaSi(i) = σ_iew^Si + κ_iew^Ca, as shown in Fig. 18 and in Table 11 for all colours. The dispersion of the full intrinsic colour can be compared to the results of previous determinations: ${\sigma }_{\mathit{BV}}^{\text{int}}=0.049$ (Jha et al. 2007), ${\sigma }_{\mathit{BV}}^{\text{int}}=0.03\pm 0.01$; ${\sigma }_{\mathit{VR}}^{\text{int}}=0.04\pm 0.01$, ${\sigma }_{\mathit{RI}}^{\text{int}}=0.06\pm 0.01$ (Nobili & Goobar 2008), ${\sigma }_{\mathit{BV}}^{\text{int}}=0.065\pm 0.008$ (Mandel et al. 2017), ${\sigma }_{\mathit{BV}}^{\text{int1}}=0.076\pm 0.011$, ${\sigma }_{\mathit{BV}}^{\text{int2}}=0.043\pm 0.013$ (2populations/decline rate, Wojtak et al. 2023). The values of ${\sigma }_{\mathit{BV}}^{\text{int}}$ found in this work for the two populations of the H_α subsample (galaxies with high and low H_α signal, as discussed in item 10 of Sect. 10) are: ${\sigma }_{\mathit{BV}}^{\text{int1}}=0.0127$, ${\sigma }_{\mathit{BV}}^{\text{int2}}=0.0149$.

Fig. 18. Correlations between the combined calcium–silicon correction DCaSi and the SALT2 variable x1 for four colours.

The smallest scatter of DCaSi is found for colour 4 ≡ B − V, with σ = 0.00349. This provides an estimate of the contribution of measurement errors to colours, all colours being processed in the same way. Extra physical variabilities worsen the evaluation of the DCaSi contribution of other colours.

10. Conclusions

The focus of this work is the extraction for each SN of our sample of extinction and intrinsic colour components, using the extinction formula for dust in our galaxy as leverage. This colour analysis is well adapted to the determination of the reddening, but does not make use of the full power of the spectral information for the intrinsic part. The linearity that we find in the correlations between different colours is nevertheless an interesting feature by itself.

1. The model presented in this work assumes the factorisation of the extinction by the host galaxy, without any correlation with the intrinsic properties of the SN. We acknowledge that the presence of circumstellar dust might invalidate this assumption.

2. The modelling of colour fluctuations requires the inclusion of four corrections to the extinction formula in the BVRI bandpasses, three intrinsic coefficients correlating the different colours, and one value for R_V. An additional scale parameter for the extinction should be set to unity. The colour accuracy achieved varies from 0.011 to 0.028 mag depending on the colour. All colours discussed in this work are ‘rest-frame’ colours. Such corrections could arise from different average dust contents in host galaxies and ours. One of the assumptions of the model (factorisation of the extinction, linearity of the correlations of intrinsic colours) could be wrong, though plausible, but we are surprised by the large difference found in the corrections to the R and I bandpasses Dex_R and Dex_I shown in Table 4 and observed directly in Fig. 5. The simulation confirms that the algorithm we use is indeed sensitive to these small corrections. The errors arising from the linear approximation cannot account for the size of the effects found, and additional intrinsic contributions are extremely unlikely to induce a correlation with extinction as strong as the one observed in Fig. 5.

3. A single extra-intrinsic contribution beyond Si IIλ4131 Å and Ca II H&Kλ3945 Å is needed to account for the colour correlations within the accuracy described immediately above We can then extract intrinsic and extinction colour components for each SN.

4. We show that the remaining magnitude fluctuation of 0.13 mag is independent of the five UBVRI bandpasses to an accuracy of better than 0.02 mag, reaching 0.007 mag for the I bandpass. There is no attempt to minimise this magnitude dispersion directly in our procedure.

5. When the extinction formula of FM99 is used as leverage to extract this extra intrinsic (one-dimensional) component, and the associated extinction component, the intrinsic coupling coefficients between the different colours are consistent over the full sample of SNe Ia at the level of 0.01–0.05. The values of the three independent coefficients are (for R_V = 2.25)

   $$\begin{array}{cc}& {\gamma }_{10}=1.265\pm 0.075,\hfill \\ \hfill & {\gamma }_{20}=1.116\pm 0.061,\hfill \\ & {\gamma }_{30}=0.539\pm 0.041.\hfill \end{array}$$(39)

   The errors quoted include the systematic errors from the colour choice, and the statistical error derived from the simulation. The extra ‘intrinsic’ content of U − I is half the size of U − B or U − V, with a large uncertainty.

6. After correcting for the analysis bias found in the simulation, the optimal mean value of the FM99 parameter R_V over our sample is found to be

   $$\begin{array}{c}\hfill R_V=2.218\pm 0.117(\text{stat.})\pm 0.050(\text{syst.}).\end{array}$$(40)

   This value is obtained within our specific treatment of the extinction scale, which is not universally adopted. The usual ‘operational’ definition R_V = A_V/E(B − V) = R_V4 then leads to R_V4 = 2.7 according to Table 4, which includes the effect of our different filter definition.

7. The associated scaling factor for the extinction component is measured to be 1.001 ± 0.044, and, as desired for consistency, is compatible with unity.

8. A previous SNfactory analysis in C11 found R_V = 2.8 ± 0.30, with a different method and a substantially different data set. A few improvements are brought here, namely the rescaling of the extinction coefficients, ensuring consistency with the B and V SN spectra; the presence of an ‘extra-intrinsic’ colour component; and we avoid an arbitrary extinction parameter for each SN (A_V). The presence of the grey fluctuation was ignored by C11, and the extraction of the extinction colour component was not explicit. The R_V value that we obtain is dependent on the method used, with an emphasis on colour correlations rather than a minimisation of Hubble residuals. The smaller number of parameters involved and the accuracy of the colour description is an argument in favour of the present method.

9. Within the accuracy of the present observations, we see no hint of a variation of the parameter R_V from one SN to another. As we find corrections to the extinction formula to be necessary, for all values of R_V in the range 1.95–2.60, the meaning of this parameter is slightly blurred. Simulations have been performed with samples including a Gaussian distribution of the value of R_V around the value of 2.20 and with an RMS of 0.5. The results do not show any significant change, and we would thus be insensitive to such a fluctuation.

10. Many instances have been given in the past of a correlation between the star formation rate (Hα signal, Rigault et al. 2013) or the host galaxy mass (Hamuy et al. 2000; Kelly et al. 2010) and the absolute magnitudes of SNe, as well as their intrinsic properties, suggesting two populations depending on their age. We checked the impact on our analysis of including a sample of 85 SNe from the SNfactory collaboration common to the Rigault et al. (2013) sample. A striking difference is seen in the amount of extinction of the two groups. The young SNe (high Hα signal of the host galaxy) have significantly more dust, with ⟨X₃⟩ = 0.0124 ± 0.036, whilst the older ones (low galactic Hα signal) have less extinction and ⟨X₃⟩= − 0.0577 ± 0.025. Our corrections slightly reduce the discrepancy in absolute magnitude of the two samples from 0.095 to 0.087 mag. The statistics of the two groups (47 and 37 SNe) are too small to evaluate a difference in the associated R_V.