© ESO 2018

1 Introduction

Since the discovery of the accelerated expansion of the Universe (Riess et al. 1998; Perlmutter et al. 1999) the world sample of Type Ia supernovae (SNe Ia) at cosmological distances has grown substantially. Mainly thanks to the combined effort from the Sloan Digital Sky Survey II (SDSS-II) and SNLS surveys, the current optical SN Ia Hubble diagram is built from 740 well measured objects (Betoule et al. 2014, and references therein) and 1049 SNe with the addition of PanSTARSS1 (Scolnic et al. 2017, and references therein). These measurements provide the most accurate study of the expansion history of the Universe to date, especially when combined with baryon acoustic oscillation (BAO) observations from 6dFGS, SDSS, BOSS and WiggleZ (Beutler et al. 2011; Anderson et al. 2014; Kazin et al. 2014; Ross et al. 2015; Alam et al. 2017) and cosmis microwave background (CMB) constraints (Planck Collaboration XIII 2016). However, it remains unclear at this point whether the precision of the measurements of cosmological parameters from, for example, LSST (Ivezic et al. 2008), WFIRST (Spergel et al. 2013), or the proposed DESIRE survey exploiting the Euclid satellite (Astier et al. 2014), most notably the equation of state parameter of the dark energy, will improve significantly with further increase in the sample of SN Ia at rest-frame optical wavelengths. Systematic uncertainties – instrumental, astrophysical and cosmological – may pose severe limitations (see, e.g., Goobar & Leibundgut 2011; Conley et al. 2011; Betoule et al. 2014; Scolnic et al. 2017).

Recent SN Ia cosmological analyses assume the same R_V for all SNe and recover significantly lower value compared to the typical Milky Way R_V ≃ 3.1. We lack solid understanding of what causes this. Furthermore, the assumption of universal R_V is almost certainly not correct. First, in the Milky Way R_V varies between 2 and 5.5 and there is no reason why this should not also be the case for other galaxies. Second, there is growing evidence that R_V varies between different supernovae (SNe) in the nearby Universe (Folatelli et al. 2010; Burns et al. 2014; Amanullah et al. 2015, and references therein). There are also indications that the most heavily reddened SNe tend to have low R_V ~ 1–2, while in the more moderately reddened events a range of values up to R_V ~ 3 is observed (Burns et al. 2014). This raises concerns for possible systematic uncertainties related to the extinction corrections. For example, that the current assumption of a universal reddening law independent of redshift may not hold. Even if the reddening law were to remain unchanged, it is already problematic that different populations would be observed as the redshift increases. At high redshift, the most reddened SNe will be missed and the sample will be dominated by low-reddening events. Applying color cuts to the nearby sample will not completely prevent the distribution of colors to differ between the low- and high-redshift sample, which can introduce systematic differences. Similarly, the “luminosity – light curve (LC) shape” relation used to calibrate SNe Ia lacks strong theoretical explanation and its evolution with redshift cannot be excluded at this stage.

One alternative that has been proposed is to explore the use of near-infrared (NIR; 1–2.5 μm) observations of SNe Ia. This has recently been used to place constraints on the Hubble constant by Dhawan et al. (2018). Observing SN Ia in the NIR has several important advantages. Firstly, the extinction corrections are by a factor of four to six smaller than in the optical B band. Secondly, there is evidence, both from observations (Meikle 2000; Krisciunas et al. 2004a,c; Wood-Vasey et al. 2008) and theoretical calculations (Kasen 2006), that at maximum SNe Ia are natural standard candles in the NIR and no correction for the LC shape is needed for normal SNe. Although more recent works by Folatelli et al. (2010) and Kattner et al. (2012) do suggest that such a correction may also be needed in the NIR, it is much smaller than in the optical. Thirdly, Krisciunas et al. (2004c,b) and later Wood-Vasey et al. (2008) demonstrated that smooth LC templates, the time-axis of which is stretched according to the SN optical stretch parameter s_opt (Perlmutter et al. 1997), accurately describe rest-frame J and H-band LCs within ten days from B-band maximum. This raises the possibility of accurately measuring the peak NIR magnitudes of SNe Ia even with a single observation within a week from B-band maximum, provided the time of B-band maximum and s_opt are known, and Barone-Nugent et al. (2012) and Weyant et al. (2014) have shown indications that this may indeed be a viable approach. If further verified, this could be very important because, due to the large atmospheric background, ground-based observations in the NIR rapidly become prohibitively time-consuming when observing faint objects and building well-sampled LCs of high-redshift SNe Ia becomes nearly impossible. The main goal of this work is to further test the feasibility to measure the NIR peak magnitudes with single observation and extend the NIR Hubble diagram to intermediate redshifts. We report on the optical UBV RI and NIR J and H measurements of 16 SNe Ia in the redshift range z = 0.037 − 0.183. We present the results of the combined analysis of this new data-set together with 102 NIR and 458 optical LCs of normal SNe Ia from the literature.

The paper is organized as follows. In Sect. 2 we present the target selection, observations, data reduction, photometry and calibration. Section 3.1 provides details on the LC fitting, including the derivation of new NIR LC templates and NIR K-corrections. In Sect. 4 we present our results and in Sect. 5 a discussion followed by our conclusions. Throughout the paper we have assumed the concordance cosmological model with Ω _M = 0.27, Ω _Λ = 0.73 and h = 0.708.

2 Observations and data reduction

2.1 Target selection and observing strategy

Between 2007 and 2008 we observed sixteen SNe Ia at redshifts between z = 0.037 and z = 0.183. Details of the targets are given in Table 1. The targets were selected from the SNe discovered by SDSS-II – 3 SNe, SNfactory (Aldering et al. 2002) – 12 SNe, and ROTSE – 1 SN, to be spectroscopically classified before maximum light. The early classification was essential in order to obtain the NIR observations within a week from maximum and to obtain good enough optical observation to derive the LC parameters. For each SN, a single NIR observation and multiple UBV RI observationswith cadence of between five and seven days were obtained, except for SNe 2007hz, 2007ih and 2007ik for which the gri photometry obtained by SDSS-II was used (Sako et al. 2018). For most SNe the observations started earlier than ten days after the B-band maximum light.

Most of the SNe in our sample are close to their host galaxy nuclei and the underlying nonuniform galaxy background may significantly degrade the photometry. To remove the galaxy background contamination we used images of the host galaxies of the SNe that we obtained about a year after the main survey, when the SN had faded. In the NIR, template images were obtained for all three instruments used. In the optical, template images wereonly obtained with ALFOSC at NOT. The galaxy template images were aligned with the SN image¹, convolved with a suitable kernel so that the point-spread functions (PSF) of the two images are the same, then scaled to match the flux level of the SN image and subtracted. The SN flux can then be correctly measured on the background-subtracted image. The image subtraction was carried out with Alard’s optimal image subtraction software (Alard & Lupton 1998; Alard 2000), slightly modified and kindly made available to us by B. Schmidt. SNe 2007hz, SNF20080512-008 and SNF20080516-022 were found to be away of their host galaxies and the galaxy background in the NIR was not a problem. For these SNe no galaxy reference images were obtained and the photometry was performed without image subtraction.

2.2 Near infrared imaging

The NIR photometry of the SNe was obtained at three facilities. J and H imaging of the three SDSS-II SNe was obtained at VLT/ISAAC (Moorwood et al. 1998) and nine of the others were observed with NOTCAM (Abbott et al. 2000) at the Nordic Optical Telescope (NOT). Six SNe were observed at the Japanese InfraRed Survey Facility (IRSF) 1.4m telescope located at SAAO and equipped with the imaging camera SIRIUS (Nagayama et al. 2003), which can simultaneously observe in JHK_s bands. Two of the SNe were observed both at IRSF and NOT. All three instruments are equipped with 1024 × 1024 pixel Hawaii detectors and filters from the Mauna Kea Observatories Near-Infrared Filter Set (MKO)² (Tokunaga et al. 2002).

The observations consisted of series of dithered images in order to facilitate the night sky background subtraction. Twilight sky flat-field images were also obtained. The image reduction followed with small modifications the scheme outlined in Stanishev et al. (2009). All ISAAC images were first corrected for the effect of “electrical ghost” with the recipe ghost within the ESO software library Eclipse.³. Then dark current and flat-field corrections were applied to all images. The sky background subtraction was performed with the XDIMSUM package in IRAF⁴. To estimate the sky background for given image, the eight exposures closest in time were scaled to the same mode and combined with median. The combined sky image was scaled to the mode of the image and subtracted from it. The pedestal on NIR arrays is highly variable in time and may be very difficult to remove accurately. In some of the images, significant residuals were found along the bottom and the middle rows of the array, where the readout of the two halves of the array starts. To correct for this effect the median value of each row was computed with 2σ clipping and subtracted. The sky-subtracted images were combined and the combined image was used to create an object mask with the IRAF OBJMASK task. The mask was then applied to each individual image and the sky subtraction was repeated, but now all pixels belonging to objects were excluded during the background estimation and the residual pedestal subtraction. The only difference was that instead of the median, the mode along the image rows was used to remove the residual pedestal.

Before the final image combination residual cosmic ray hits were identified and cleaned with the Laplacian detection algorithm ofvan Dokkum (2001). ISAAC and NOTCam have significant optical image distortions. We derived corrections for this effect using observations of dense stellar fields with astrometry from SDSS and combined the images with the DRIZZLE algorithm (Fruchter & Hook 2002). For the final combination, the individual images were also assigned weights in order to optimize the signal-to-noise ratio (S/N) for faint point sources. The weights are inversely proportional to the product of the square of the seeing, the variance of the sky background and the multiplicative factor that brings the images to the same flux scale.

Photometry and calibration. The instrumental magnitudes were measured by PSF photometry with the DAOPHOT package (Stetson 1987) in IRAF. The photometric calibration of the SNe observed at ISAAC was obtained with standard stars from Persson et al. (1998), which were observed as part of ISAAC standard calibration plan. The zero-points (ZPs) estimated from the standards were transferred to the SN images assuming atmospheric extinction coefficients of k_J = 0.05 and k_H =0.04 mag airmass⁻¹ (Lombardi et al. 2011). The image ZPs for the observations obtained at NOT and SAAO were computed with 2MASS stars in the field of view. For each SN image, a ZP was calculated for each 2MASS star. The image ZPs and their uncertainty are, respectively, the average of the individual ZPs (with 3σ outliers removed if present) and the standard deviation of the mean. No color terms were applied for the following reasons. Firstly, it is quite difficult to derive accurate NIR color-terms observationally because of the lack of suitable standard star fields containing stars with different colors, for example, like those used in the optical (Landolt 2009). Secondly, for the photometricsystem used we expect small color terms. Synthetic photometry of the stars in Pickles (1998) with the passbands that we adopted for the MKO and 2MASS systems (see the Appendix) shows only small color-terms (≤0.03) between these two systems. Lombardi et al. (2011) also found zero color-term for ISAAC H band and only a small one ~0.02 (J − H) for J band. The final calibrated magnitudes are given in Table 2.

2.3 Optical imaging

For SNe 2007hz, 2007ih and 2007ik the photometry obtained by SDSS-II was used (Sako et al. 2018). For the remaining 13 SNe, the optical observations were obtained with four different instruments, ALFOSC, MOSCA and STANCam at NOT equipped with broadband UBV RI filters and ANDICAM at CTIO 1.3m telescope operated by the SMARTS Consortium equipped with BV RI filters. Standard fields from Landolt (2009) were also observed with each instrument. The observations were collected between May and July 2008. On May 23rd, 2009 the SN fields were re-observed with ALFOSC to obtain reference images for the galaxy subtraction. All CCD images were bias and flat-field corrected, and trimmed. Cosmic ray hits were identified and cleaned with the Laplacian detection algorithm of van Dokkum (2001). When multiple exposures per filter were obtained, they were combined into a single image.

Photometry and calibration. For the optical observations, aperture photometry was chosen over PSF photometry because the main instrument used for the survey – ALFOSC – has significant PSF variation over the image and in most SN fields there were not enough bright stars to model this effect. The instrumental magnitudes were measured through digital apertures with diameter of double the mean full width at half maximum (FWHM) of the stars in the field. It is known that aperture diameter of approximately FWHM maximizes the S/N, but double the FWHM aperture was chosen in order to minimize the effect of PSF variation.

The standard star fields observed in eight photometric nights with ALFOSC at NOT were used to calibrate local sequences of stars around the SNe. Aperture corrections from two- to five- times FWHM were computed with the IRAF.PHOTCAL package and applied to the measured magnitudes. All standard star measurements were fitted simultaneously (with 3σ clipping) with linear equations of the form (Harris et al. 1981): $$\begin{array}{lll}U\hfill & =\hfill & u+c{t}_U(U-B)-{k^{\prime }}_{U}X-{k^{″}}_{U}X(U-B)+z{p}_U\hfill \\ B\hfill & =\hfill & b+c{t}_B(B-V)-{k^{\prime }}_{B}X-{k^{″}}_{B}X(B-V)+z{p}_B\hfill \\ V\hfill & =\hfill & v+c{t}_V(B-V)-{k^{\prime }}_{V}X+z{p}_V\hfill \\ R\hfill & =\hfill & r+c{t}_R(V-R)-{k^{\prime }}_{R}X+z{p}_R\hfill \\ I\hfill & =\hfill & i+c{t}_I(V-I)-{k^{\prime }}_{I}X+z{p}_I,\hfill \end{array}$$(1)

where upper-case and lower-case letters denote the standard and instrumental magnitudes, respectively, ct the color-terms, X the airmass, k′ and k^″ the first- and the second-order extinction coefficients and the zp the instrument ZPs. The derived instrument color-terms are given in Table 3. The second-order extinction coefficients k^″ were fixed to −0.03 based on our synthetic photometry calculations and past photometric works (see also Bond 2005). Because the observations were obtained during a short period of time the data was fitted with common color-terms and instrument ZPs. The standard star observations with the other three instruments we only used to derive the color-terms, again assuming k^″ = −0.03 for U and B bands.

With the derived transformation equations, the standard magnitudes of the stars in the SN fields were computed. The weighted mean and standard deviation of the calibrations from different nights (weighted by the estimated photometricuncertainties) were taken as standard magnitude and its uncertainty. The stars that showed scatter larger than that expected from the photometric errors were discarded. The final calibrated local sequences are indicated onthe finding charts in Fig. A.1 and the magnitudes are listed in Table A.1.

The calibrated local sequence stars, measured with two-times FWHM diameter apertures, were used to compute the image ZPs. For each SN image, a ZP zp^X ≡ zp − k′ X was calculated for each calibrated star by applying Eqs. (1). The final image ZPs and their uncertainty are, respectively, the average of the individual ZPs (with 3σ outliers removed if present) and the standard deviation. The image ZPs were added to the SN magnitudes measured in the galaxy-subtracted images to obtain the SN magnitudes in the natural systems of the instruments used, m^nat.

It is known that Eqs. (1) cannot accurately transform SN photometry into the standard Johnson-Cousins UBV RI photometricsystem (Suntzeff 2000; Stritzinger et al. 2002; Krisciunas et al. 2003). For this reason, the optical photometry was calibrated with the method described in Stanishev et al. (2007), which is based on the S-correction introduced by Stritzinger et al. (2002). Because we took into account the second-order extinction, the calibration equations were slightly modified. The calibration correction from natural magnitudes at airmass X, $m_X^{\text{nat}}$, to standard magnitudes m^std is computed by means of synthetic photometry: $$\begin{array}{lll}{m}^{\text{std}}\hfill & =\hfill & m^{\text{nat}}-2.5\log \left({\displaystyle \int }{f}_{\lambda }^{\text{phot}}(\lambda )R^{\text{std}}(\lambda )\text{d}\lambda \right)\hfill \\ \hfill & \hfill & +2.5\log \left({\displaystyle \int }{f}_{\lambda }^{\text{phot}}(\lambda )R^{\text{nat}}(\lambda )p{(\lambda )}^X\text{\hspace{0.17em}d}\lambda \right)+\text{const}.\hfill \end{array}$$(2)

Here $f_{\lambda }^{\text{phot}}(\lambda )$ is the photon spectral energy distributions (SED) of the SN, R^nat(λ) and R^std(λ) are the responses of the natural and the standard systems, respectively. We note that R^nat(λ) does not include the atmospheric transmission and thus the transmission of the atmosphere at airmass X = 1 above the site of observation, p(λ), is also included in the calculation. Because the SNe in our sample do not have extensive spectrophotometry we used Hsiao et al. (2007) spectral template as $f_{\lambda }^{\text{phot}}(\lambda )$, corrected for the Milky Way reddening and red-shifted to the redshift of each SN. The constant in Eq. (2) is such that the correction is zero for A0 V stars with all color indexes zero. This ensures that for normal stars the synthetic S-correction gives the same results as the linear color-term corrections (Eq. (1)). The constant can be determined from synthetic photometry of stars for which both photometry and spectrophotometry is available. In the Appendix we provide full account of the calibration procedure, including a discussion of the uncertainties introduced when the second-order extinction is not taken into account. The final photometry of the SNe is given in Table 4.

3 SNe Ia as standardizable candles

The rest-frame peak magnitudes SN Ia in band X, m_X, can by standardized with the following two-parameter relation: $$m_X^{\ast }=m_X-a_X\times (\Delta M_{15}-1.1)-R_X^{YZ}\times E(Y-Z),$$(3)

where E(Y − Z) ≥ 0 is the color excess in bands Y and Z caused by dust and which can be converted into absorption in band X via $R_X^{YZ}\times E(Y-Z)$. The second correction term in the equation accounts for the well-known “luminosity – LC shape” relation and as a LC shape parameter we use Δ M₁₅, which measures the rest-frame B-band magnitude decline during the first 15 days after the time of B-band maximum light $t_{B_{\text{max}}}$ (Phillips 1993). Throughout this work we use E(B − V), which is calculated from the observed color excess at the time of B-band maximum light $E{(B-V)}_{\text{obs}}={(B-V)}_{B_{\text{max}}}-{(B-V)}_{B_{\text{max}},0}$ by setting the negative values to zero. The intrinsic unreddened SN Ia color ${(B-V)}_{B_{\text{max}},0}$ is a nonlinear function of ΔM₁₅ that we determined (see Appendix F).

3.1 Measuring the light curve parameters

SN Ia observations rarely allow direct measurement of the peak magnitudes and Δ M₁₅. Thus, to obtain these parameters a multitude of indirect methods that fit optical and/or NIR template LCs to the observations have been developed (Perlmutter et al. 1997; Goldhaber et al. 2001; Guy et al. 2005, 2007; Jha et al. 2007; Conley et al. 2008; Prieto et al. 2006; Krisciunas et al. 2004b; Wood-Vasey et al. 2008; Mandel et al. 2009, 2011; Burns et al. 2011). For the needs of our analysis we developed our own stretch parametrization-based fitter (Perlmutter et al. 1997; Goldhaber et al. 2001) for the following reasons. First, our primary goal was to test the accuracy of the estimate of SN Ia NIR peak magnitudes from a single observation obtained close to maximum. Krisciunas et al. (2004b) have presented evidence that the stretch approach could be applied to the first maximum of J band LCs and that the NIR and the optical stretch are tightly correlated. With the stretch fitting approach we were able to verify this result with a much larger sample of SNe with well-sampled LCs. Second, we used our own NIR templates and prescription to compute NIR K-corrections (Appendix E and 3.1.2). Third, from the optical LCs we need $t_{B_{\text{max}}}$, B and V peak magnitudes, and a LC shape parameter. Considering the above, stretch-based fitter appears to be best suited for our purposes.

Our fitting program is based on the IDL fitting package MPFIT⁵. It can fit multicolor LCs simultaneously using two approaches. The first is to fit the data with rest-frame template LCs coupled with a spectral template to compute K-corrections (Kim et al. 1996) in order to transform from rest-frame to observed magnitudes. In the second approach a spectral template is red-shifted and directly integrated through the observed filters. This requires that the spectral template is calibrated to match the rest-frame LCs of SNe Ia, which is the case for Hsiao et al. (2007) template. This template is generated in such a way that it matches the UBV RI LC templates of Knop et al. (2003) and its B − V color at maximum is − 0.058 mag. In both approaches the fitted parameters are the common stretch s that multiplies the time axis of the template, the time of B-band maximum light $t_{B_{\text{max}}}$, and the magnitude offsets needed to fit each SN LC. The full covariance matrix of the fit parameters is also computed. The LC templates are normalized to zero magnitude at $t_{B_{\text{max}}}$ and becauseeach template magnitude offset is fitted independently, the procedure returns the rest-frame SN magnitudes at $t_{B_{\text{max}}}$. The program is flexible and one or several fit parameters can be fixed. It also allows the spectral template to be reddened by arbitrary reddening laws.

The stretch parametrization can only be used for U, B or V bands. In the redder bands SNe Ia show a second maximum (or just a shoulder in R) and the stretch approach can only be applied to the first maximum (e.g., Krisciunas et al. 2004b). Therefore, the fits were limited to the observations obtained up to phase + 10 stretch-corrected days in J and H bands, and up to phase + 35 days in the optical. The stretch-corrected SN phase, ϕ, is defined as $$\varphi =C\frac{t-t_{B_{\text{max}}}}{s(1+z)},$$(4)

where z is the SN host galaxy redshift, s is the stretch and the constant C is set to 1, unless specified otherwise. The stretch parameter s can be either the optical stretch s_opt or the NIR stretch s_NIR, which will be defined later.

3.1.1 Optical light curve parameters

For the analysis of the NIR LCs we need the optical stretch s_opt, $t_{B_{\text{max}}}$, and B and V peak magnitudes. To determine them, the Milky Way extinction-corrected g and r LCs of SNe 2007hz, 2007ih and 2007ik by SDSS-II, and the B and V -band LCs of the remaining objects werefitted simultaneously⁶ in a two-iteration procedure. In the first iteration the Hsiao et al. (2007) spectral template was integrated through the relevant filters and fitted to the data to obtain the SN B − V color at $t_{B_{\text{max}}}$, ${(B-V)}_{B_{\text{max}}}$. The spectral template was then modified to match the observed SN ${(B-V)}_{B_{\text{max}}}$ color⁷ and the fitting was repeated to obtain the final parameters. This two-step procedure has the advantage that it accounts for the LC shape modification by the redshift of the filters’ effective wavelength caused by reddening. In particular, the derived stretch parameters are largely free from this effect. The optical LC parameters of the SNe from our new sample are shown in Table 7. They show that there are no underluminous or highly reddened SNe in our sample. Further, the optical observations of all but two SNe started before phase +15 days and the NIR observationsof only two SNe were obtained after phase +10 days. For illustration purposes, the simultaneous gri or UBV RI fits are shown in Fig. 1. In the analysis we use the results from the simultaneous fits of the g and r, and B and V bands only.

The optical LCs for the literature NIR sample (102 objects) were taken from the following sources: Calàn-Tololo (Hamuy et al. 1996), CfA1,2,3 and 4 (Riess et al. 1999; Jha et al. 2006; Hicken et al. 2009a, 2012), Berkeley (Ganeshalingam et al. 2010), Carnegie Supernova Project (CSP; Contreras et al. 2010; Stritzinger et al. 2011), Supernova Cosmology Project (SCP; Kowalski et al. 2008), K. Krisciunas and collaborators (Krisciunas et al. 2004c,b, 2003, 2001, 2000, 2007), SN 1998bu (Hernandez et al. 2000; Jha et al. 1999; Suntzeff et al. 1999), SN 2000E (Valentini et al. 2003), SN 2002bo (Benetti et al. 2004), SN 2002dj (Pignata et al. 2008), SN 2003cg (Elias-Rosa et al. 2006), SN 2003du (Stanishev et al. 2007), SN 2005cf (Wang et al. 2009b; Pastorello et al. 2007b), SN 2008Q (Stanishev et al., in prep.), and SN 2011fe (Munari et al. 2013). For the SNe observed in the NIR by Barone-Nugent et al. (2012) (hereafter referred as BN12 sample) the optical g and r LCs from Maguire et al. (2012) were used. The CSP photometry is given in the natural photometric systems of the CSP telescopes and the observed B and V bands used were those from Stritzinger et al. (2011). For the other SNe the standard Bessell B and V bands (Bessell & Murphy 2012) were used. If a given SN had optical LCs available from different sources, the optical parameters associated to it were derived from the LCs that had better phase coverage and sampling. In a few cases it was necessary to combine the optical photometry in order to obtain more reliable LC fits.

To put the optical properties of the SN in our new sample in the context of other SNe Ia, all objects with LCs that allow stretch fitting were also analyzed. This analysis was also performed because it may be relevant and of practical interest for SN cosmology. For many current and future high-redshift SN samples only the rest-frame B and V band are in the optical window. The redder bands are red-shifted in the NIR and are or will be unavailable. The UV part of the spectrum is red-shifted to the optical, but the UV properties of SNe Ia are much less understood and the peak magnitudes appear to have larger intrinsic scatter (Brown et al. 2010; Milne et al. 2010, 2013). Given this, it is interesting to study how well the SN Ia standard candle can be calibrated when only rest-frame B and V observations are available, now with much larger sample of nearby SNe. In addition, with this large sample of high-quality LCs we can study the SN intrinsic B − V color at $t_{B_{\text{max}}}$, ${(B-V)}_{B_{\text{max},0}}$, which is important for estimating the dust extinction in the SN host galaxies.

In our analysis we found that the “luminosity – LC shape” relation for normal SNe Ia is linear (within the scatter) when Δ M₁₅ is used as a LC shape parameter and not the stretch. Using a subsample of ~200 SNe with well-sampled LCs we derived a new s_opt −ΔM₁₅ relation for normal SNe Ia (see Appendix D for details) and used it to convert stretch to Δ M₁₅.

In Table 5 we report the weighted mean difference and 1σ scatter between the B and V peak magnitudes, ${(B-V)}_{B_{\text{max}}}$, and Δ M₁₅ for the data sets that have more than ten SNe in common. One can see that there are small systematic differences in the peak magnitudes, which are significant to about 1σ level. Similar small differences between the diferent datasets have also been noticed previously, for example, Ganeshalingam et al. (2013). The reasons for these systematic differences are likely related to the photometric calibration and different LC sampling. ${(B-V)}_{B_{\text{max}}}$ and Δ M₁₅ shows relatively large scatter, but no significant bias. We also compare our results those obtained by Betoule et al. (2014) with the SALT2 LC fitter. For the 104 SNe in common the SALT2 x1 and our stretch parameters are tightly related as s_opt = 1.0 + 0.108 × x1, and B_max show only a small offset with a median value of 0.022 mag and a robust 1σ scatter of 0.042 mag. The median difference between our ${(B-V)}_{B_{\text{max}}}$ color and SALT2 c parameter is 0.047 ± 0.042 mag, which is consistent with the definition of the c parameter (Guy et al. 2005). Given the significant differences between the two LC fitting approaches, we conclude that the results compare rather well.

Fig. 1 griJH and UBV RIJH LCs and the best fits. For the three SDSS-II SNe (2007hz, 2007ih, 2007ik) from bottom to top are plotted g, r − 0.4, i − 1, J − 2.8 and H − 3.8 (no J for SN 2007ik). For the other 13 SNe are plotted U + 0.8, B, V − 0.2, R − 0.6, I − 1.6, J − 3 and H − 3.5. We note that the simultaneous gri or UBV RI fits were only used for these plots. In the analysis we use the results from the simultaneous fits of the g and r, and B and V bands only.

Fig. 2 Data used to create the J and H LC templates. The templates areoverplotted with solid lines. σw is the weighted rms of the data around the template during the phases earlier than +10 days.

3.1.2 NIR light curve parameters

The literature NIR LCs used in this work come from the following sources: BN12, Wood-Vasey et al. (2008) (CfA set), Contreras et al. (2010); Stritzinger et al. (2011) (CSP set), Krisciunas et al. (2000, 2001, 2003, 2004c,b, 2007) (Kri set) and the following SNe from various sources (Var set) SN 1998bu (Hernandez et al. 2000; Mayya et al. 1998; Jha et al. 1999), SN 1999ee (Hamuy et al. 2002), SN 2000E (Valentini et al. 2003), SN 2002dj (Pignata et al. 2008), SN 2003cg (Elias-Rosa et al. 2006), SN 2003du (Stanishev et al. 2007), SN 2003hv (Leloudas et al. 2009), SN 2004eo (Pastorello et al. 2007a), SN 2005cf (Wang et al. 2009b; Pastorello et al. 2007b), SN 2008Q (Stanishev et al., in prep.), and SN 2011fe (Matheson et al. 2012).

NIR light curve templates

For the fit of the NIR LCs we adopted the template fitting approach coupled with a procedure to compute NIR K-corrections that we developed (see Appendix E for details). The K-corrections were computed from the different observed filter systems to the Mauna Kea Observatories JHKs photometricsystem (Tokunaga et al. 2002) as realized by the filter - detector combination at the NOTCam instrument at the NOT. Although there are currently several NIR LC templates available, with the excellent new NIR photometry of many nearby SNe Ia published recently (e.g., Contreras et al. 2010; Stritzinger et al. 2011) we also derived new templates.

To derive the NIR template LCs, only SNe that had (i) good optical LCs and (ii) several NIR observations obtained before phase +10 days and well-distributed along the LC were used. Known peculiar objects such as SN 2002cx-likes, SNe 2006bt, 2006ot, and SN 1991bg-likes, were discarded. As a starting point, the J and H-band LCs were simultaneously fitted with the Wood-Vasey et al. (2008) templates and our new NIR K-correction procedure. In the fit, the templates were stretched with the optical stretch s_opt and the zero time was fixed to $t_{B_{\text{max}}}$ as in Krisciunas et al. (2004c,b). Examination of the well-sampled LCs revealed that some fits were not very good. The reason appeared to be that the optical stretch s_opt did not provide adequate stretching to fit the NIR LCs. To further test this we selected a subsample of SNe with well-sampled LCs, which we refer to as the “GOLD” sample. For a SN to enter this sample it has to have (i) at least four points in J and/or H band obtained before phase +10 days with small uncertainty (in general < 0.05 mag); (ii) the points are evenly-distributed along the LC; (iii) has observations within two days from the time of NIR maximum and at least one observation between phases 5–10 days. Such SNe allow to perform the fit of the first NIR LC maximum with the NIR stretch s_NIR and the time of maximum as free fitting parameters. The 40 objects in the GOLD sample are: SN 1998bu, SN 1999ee, SN 1999ek, SN 2000E, SN 2001cz, SN 2001el, SN 2002bo, SN 2002dj, SN 2003du, SN 2004eo, SN 2004ey, SN 2005M, SN 2005cf, SN 2005el, SN 2005eq, SN 2005hc, SN 2005iq, SN 2005kc, SN 2005ki, SN 2006D, SN 2006X, SN 2006ax, SN 2006et, SN 2006kf, SN 2006le, SN 2006lf, SN 2006os, SN 2007A, SN 2007S, SN 2007af, SN 2007bc, SN 2007bm, SN 2007ca, SN 2007le, SN 2008Q, SN 2008bc, SN 2008fp, SN 2008gp, SN 2008hv, and SN 2011fe. SNe 2005el, 2005eq, 2005iq, 2006X, and 2006ax have LCs from both CSP and CfA, and the LCs were analyzed separately. The fit of the CfA LC of SN 2005iq and SN 2006ac provided very weak constraints on the stretch and the time of maximum and was moved from the GOLD sample to the SILVER one (see Sect. 3.1.3). The differences between the peak magnitudes estimated from the CSP and CfA photometry are at most 0.08 mag. The only exception is H band of 2006ax, for which the H-band photometry from CfA was 0.22 mag fainter than CSP. Even though the analysis of the maximum colors indicated that those obtained from the CSP photometry were more likely to be correct, we also kept the CfA photometry of this SN.

To derive the final templates, the observed LCs were K-corrected to rest frame magnitudes, normalized to zero magnitude at $t_{B_{\text{max}}}$, converted to phases using Eq. (4) with s = s_NIR, and stacked together. The data for phases before + 10 days were fitted with smoothing-spline function, manually adjusting the smoothing parameter to achieve a smooth fit that leaves no correlated residuals and again taking care to avoid overfitting. A few points that deviated by more than 3σ were removed from the fits. The same was repeated for the later data and the two fits were stitched together. Because our NIR K-correction are reliable up to phases ~+40–50 days, the data past + 50 days was not used. We noticed that the NIR stretch s_NIR of SNe with s_opt ~ 1 was not 1. For this reason the template time axis was additionally stretched so that the SNe with s_opt~ 1 also had on average s_NIR ~ 1.

Armed with the new templates the above steps were repeated to derive our final J and H-band LC templates shown in Fig. 2. The templates were additionally normalized to have zero magnitude at the first NIR LC peak so that the fitted magnitude offsets are the rest-frame NIR peak magnitudes. The weighted root mean square (rms) scatter of the points around templates⁸, σ_w, is σ_w ≃ 0.02 mag and σ_w ≃ 0.03 mags for J and H, respectively. The small scatter shows that the templates describe the LCs of the GOLD sample very well. Between +10 and +50 days the scatter is considerably larger, ~0.13 mag and 0.09 mag for J and H, respectively.

Fig. 3 Panel a: NIR vs. optical stretch parameter. The solid line is the best-line fit. The four encircled points are the SNe shown in Fig. 4; Panel b: histogram of the difference between time of B band maximum estimated from the NIR and optical LC fits from the GOLD sample.

Fig. 4 Example of NIR LC fit of SNe from the GOLD sample. The solid line is the best fit, and dotted lines are fits with sNIR = sopt. The J band is shown in blue and H in red.

GOLD SN sample

The LCs of the GOLD SN sample were fitted with the new templates with the NIR stretch s_NIR and the time of maximum as free fitting parameters. The NIR and the optical LC parameters are shown in Table 6. Figure 3a shows the NIR vs. the optical stretch parameters. The best weighted least-squares linear fit to the points is shown with the solid line and it essentially coincides with the one-to-one correspondence s_NIR = s_opt. An inspection of the NIR LC fits suggests that the scatter in Fig. 3a most likely arises from real differences between optical and NIR stretch. In Fig. 4 are shown the fits of four SNe. The solid lines are the fits with s_NIR as a free parameter and the dotted lines with s_NIR = s_opt. It is evident that the LCs of SN 2002bo and SN 2006et are better fitted with s_NIR as a free parameter. In order to have χ²∕DoF = 1 (DoF− degrees of freedom) with respect to the line defined by the s_NIR = s_opt relation it was necessary to add the quantity σ_s = 0.05 in quadrature to the uncertainties of s_NIR. This can be interpreted as intrinsic scatter of the s_NIR − s_opt relation.

The zero time of the NIR LC templates is set at $t_{B_{\text{max}}}$. Thus, the NIR LC fits with the stretch and time of “maximum” as free parameters provide an independent estimate of $t_{B_{\text{max}}}$. In Fig. 3b we show the histogram of the difference $\Delta t_{B_{\text{max}}}$ between the estimates of $t_{B_{\text{max}}}$ from the NIR and the optical LC fits. The distribution peaks at zero and has a FWHM around one day. The width is larger than what can be expected from the uncertainties of $\Delta t_{B_{\text{max}}}$ alone and indicates possible presence of small random time-offsets between the optical and NIR LCs. Thus, an additional term σ_t was added in quadrature to the uncertainties of $\Delta t_{B_{\text{max}}}$ and its value σ_t = 0.44 day was estimated so that χ²∕DoF with respect to the weighted mean of $\Delta t_{B_{\text{max}}}$ was equal to one.

3.1.3 Fitting under-sampled NIR light curves

The analysis of the GOLD sample further reinforced the results of Krisciunas et al. (2004c,b) and Wood-Vasey et al. (2008) who showed that smooth LC templates stretched according to the optical stretch accurately describe the rest-frame J and H-band LCs within ten days from B-band maximum. Thus, if s_opt and $t_{B_{\text{max}}}$ are estimated from optical observations, by fixing the zero time of the NIR template to $t_{B_{\text{max}}}$ and setting s_NIR = s_opt, one could inprinciple measure the SN Ia NIR peak magnitudes with a single observation. We used this approach to estimate the peak magnitudes of the SNe in our sample and all published SNe with under-sampled LCs with observationsbefore +10 days. This also includes the BN12 sample. However, we note that most SNe in this sample have severalobservations around maximum and the LC quality is closer to that of the GOLD sample; often the only difference is in the photometric uncertainties, which are smaller in the GOLD sample. The parameters derived from the LC fits are given in Tables 7–9.

Fig. 5 Peak magnitude uncertainties of a SN with sNIR = 1.0 vs. the phase of the single observation. The 1σ uncertainties of the stretch and the time of maximum were assumed to be 0.05 and 0.44 days, respectively.

Expected uncertainties

With the adopted fitting procedure the uncertainties of the time of maximum and the stretch, both from the optical LC fitting and the nonzero scatter σ_t and σ_s, will propagate to the uncertainties of the estimated peak magnitudes. This additional source of uncertainty is especially important when only a single NIR observation is available. In this case the magnitude of the uncertainty will also depend of the phase of the observation.

To estimate the expected magnitude of this uncertainty when only a single observation is available a series of simple Monte Carlo simulations were performed. For a given stretch s₀, 10 000 random pairs of stretch and time offset (s_0,i, t_0,i) were generated with the relations $s_{0,i}=s_0+N(0,{\sigma }_s^2)$ and $t_{0,i}=N(0,{\sigma }_t^2)$, where N(μ, σ²) is Gaussian distribution with mean μ and variance σ², and σ_t = 0.44 day and σ_s = 0.05. The J and H templates were stretched and shifted with each pair (s_0,i, t_0,i), and the scatter with respect to the unperturbed template was computed for each phase from − 7 to +10 days. The results for s_NIR = 1.0 are shown in Fig. 5. As expected, when the observation is obtained with a few days from the NIR maximum the peak magnitudes in both bands are fairly insensitive to the uncertainties of the stretch and the time of maximum. After maximum, the uncertainty in the H band remains less than 0.02 mag, while in J it increases reaching 0.12 mag at +10 days. These results can be understood by noting that the J-band flux changes much more rapidly than H.

The previous simulations do not include the uncertainties from the optical LC fitting. Thus, a set of more realistic simulations based on the LCs of the GOLD SN sample were also performed. Two points were randomly drawn from each LC, one from the data within 3 days from the NIR maximum, and another one from the data taken later that 4 days after $t_{B_{\text{max}}}$. The selected points were fitted with the method for under-sampled LCs, and the weighted mean and RMS of the difference from the peak magnitudes estimated from the full LC fits was computed. The mean values were in all cases close to zero indicating that the LC under-sampling does not lead to bias of the estimated peak magnitudes. The weighted scatter for the points close to maximum was ~0.025 mag in both bands. The scatter increases considerably with the points obtained after phase +4 days: 0.11 mag in J and 0.07 mag in H. Based on the results of the simulations the literature sample with under-sampled LCs was split into two groups (i) SILVER with the SNe that have observations within two to three days from the expected time of the NIR maximum light, and (ii) BRONZE with the remaining objects. Table 10 shows statistics of the NIR SNe samples.

Barone-Nugent et al. (2012) and Weyant et al. (2014) also performed simulations to estimate the uncertainties from under-sampling the NIR LCs. Although these authors used different approaches from ours, the conclusions was similar - with prior information from optical LCs, under-sampled NIR LCs can provide accurate enough peak magnitudes to be useful in cosmology.

Propagating the uncertainty

During the fitting the uncertainties of the optical LC stretch and time of B-band maximum were propagated to the uncertainty of the peak magnitudes with the following Monte Carlo procedure. The stretch and the time of maximum were perturbed according to the relations $s=s_{\text{opt}}+N(0,{\sigma }_s^2+d{s}_{\text{opt}}^2)$ and $t=t_{B_{\text{max}}}+N(0,{\sigma }_t^2+d{t}_{B_{\text{max}}}^2)$ and the LC was fitted with the perturbed pair (s, t). Here ds_opt and $d{t}_{B_{\text{max}}}$ are the uncertainties of the stretch and time of maximum, and σ_t = 0.44 days and σ_s = 0.05 as determined previously. For each SN this was repeated 10 000 times and the standard deviation of the resulting peak magnitudes was added in quadrature to the uncertainty derived by the LC fitter. The covariances between s_opt, $t_{B_{\text{max}}}$, and the peak magnitudes were also computed from the simulations.

Fig. 6 Panel a: Jmax − Hmax color vs. ΔM15. The solid line shows the least-squares linear fit to the GOLD sample; panel b: Jmax − Hmax $vs.E{(B-V)}_{\text{obs}}$. The solid line is the least-squares linear fit to the GOLD sample, the dashed and the dotted lines are the expected relation for Milky Way-type dust and RV = 1.7 and 3.1, respectively.

4 Results

4.1 J_max − H_max, $V_{B_{\text{max}}}-J_{\text{max}}$ and $V_{B_{\text{max}}}-H_{\text{max}}$ colors

Figure 6 shows the J_max − H_max color vs. ΔM₁₅ and $E{(B-V)}_{\text{obs}}$. The J_max − H_max color shows a fair amount of scatter, largely due to our sample and the BRONZE sample. However, the scatter in the GOLD sample is fairly small. The weighted least-squares linear fit to J_max − H_max vs. Δ M₁₅ in the GOLD sample only indicates a possible weak dependence on ΔM₁₅ with a slope of 0.14 ± 0.06. The fit also implies that a typical SN Ia with ΔM₁₅ = 1.1 has J_max − H_max = −0.23 ± 0.02 mag. For χ²∕d.o.f. of the fit to be one it is necessary to add 0.08 mag in quadrature to the errors, which can be interpreted as the intrinsic scatter of the J_max − H_max color. Similar intrinsic scatter is also inferred from the BN12 sample, however, J_max − H_max color in this sample is systematically redder by about 0.1 mag.

The GOLD sample also reveals a weak linear relation between J_max − H_max and $E{(B-V)}_{\text{obs}}$. The relation is shown with a solid line in Fig. 6b. The expected relations for Milky Way-type dust with R_V = 1.7 and 3.1 are shown with the dashed and the dotted lines, respectively. The slope of the best-fit line, 0.12± 0.03, is consistent with a low R_V value around 1.5. The fit indicates that an un-reddened SN should have J_max − H_max = −0.25 ± 0.03. We also note that there is no apparent dependence of J_max − H_max of the redshift. For both the GOLD sample with mean redshift z ≃ 0.02 and for our sample with mean redshift z ≃ 0.09 J_max − H_max color is − 0.24 mag.

Our estimate of the intrinsic J_max − H_max color of a SN with ΔM₁₅ = 1.1 mag is similar to the results of Krisciunas et al. (2004c) and Wood-Vasey et al. (2008). It is also consistent with Burns et al. (2014) who, using a different approach, found that a SN with stretch equal to one has J_max − H_max = −0.23 mag with intrinsic dispersion of 0.09 mag. It is interesting to note, however, that Folatelli et al. (2010) based on part of the data we used in our analysis find it to be close to zero. From Fig. 6 one can see that this is inconsistent with our results. We also note that the observed variation of the J_max − H_max color is far too big to be caused by dust extinction of any reasonable law.

Finally, Fig. 7 shows $V_{B_{\text{max}}}-J_{\text{max}}$ and $V_{B_{\text{max}}}-H_{\text{max}}$ colors vs. $E{(B-V)}_{\text{obs}}$. Most of the SNe that have $E{(B-V)}_{\text{obs}}\ge 0.4$ mag align along the line tracing Milky Way-type extinction with R_V = 1.7 − 1.9. There is also a set of SNe with $E{(B-V)}_{\text{obs}}$ between 0.1 mag and 0.6 mag that seem tofollow the line with R_V = 3.1. The two SNe with $E{(B-V)}_{\text{obs}}<-0.15$ mag appear to be outliers. However, one of them SNF20080522-004 has large uncertainty of the maximum B − V colors of ~0.2 mag. The other object is SN 2006hx, which is also an ~1 mag outlieron the optical Hubble diagram, hence its anomalous V −NIR colors.

4.2 Hubble diagram

The main goal of SNe Ia Hubble diagram analysis is to derive the cosmological parameters, such as the mass density of the Universe. This requires a sample of SNe spanning from low- to high-redshifts (Goobar & Perlmutter 1995). The nearby SNe Ia sample alone cannot constrain the cosmological parameters but instead can be used to derive the calibration parameters a_X and R_X in Eq. (3), the absolute SN Ia magnitudes and to assess the precision of the SN Ia standard candle. In our analysis the a_X and R_X parameters, and the absolute magnitude in the band X, M_X, of a fiducial SN Ia with Δ M₁₅ = 1.1 mag were determined by χ²-minimization of: $${\chi }^2={\sum }_{i=1}^{N_{\text{set}}}{\sum }_{j=1}^{N_{\text{SN}}}\frac{{\left(m_{X,i,j}^{\ast }-M_X^i-{\mu }_{i,j}\right)}^2}{{\sigma }_{\text{lc},i,j}^2+{\sigma }_{m_X,\text{intr},i}^2+w{\sigma }_{\text{ecorr},i,j}^2+{\sigma }_{\text{ext},i,j}^2}.$$(5)

with iterative 3σ outlier rejection and where the distance modulus μ was computed with the adopted cosmological parameters. In Eq. (5) we allowed for the different data sets to have different absolute magnitudes $M_X^i$. Here σ_ext account for the random host galaxy peculiar velocities with dispersion σ_pec = 150 km s⁻¹ and ${\sigma }_{m_X,\text{intr}}$ accounts for possible intrinsic SN Ia luminosity scatter, which cannot be corrected for by the model defined by Eq. (3). In principle, ${\sigma }_{m_X,\text{intr}}$ should also include contribution from any other sources of uncertainty such as calibration uncertainties, possible mismatch between the templates and the actual fitted LCs, etc., which could be different in different data sets. For this reason, the values of the terms ${\sigma }_{m_X,\text{intr}}$ were obtained so that the reduced χ² is one for each data set. This approach helps to prevent samples with poorer-quality data from decreasing the accuracy of the derived parameters (Amanullah et al. 2010). For low-redshift samples the inferred value of ${\sigma }_{m_X,\text{intr}}$ is also dependent on the adopted value of the dispersion of the galaxies peculiar velocity σ_ext, which is somehow uncertain. The term σ_lc $${\sigma }_{\text{lc}}^2={\sigma }_{m_X}^2+a_X^2{\sigma }_{\Delta M_{15}}^2-2a_X{\sigma }_{[\Delta M_{15},m_X]}.$$(6)

accounts for the uncertainties of the peak magnitude and Δ M₁₅ and their covariance derived from the LCfit. Hereafter σ_[X,Y] denotes co-variance of the variables X and Y. The term σ_ecorr $$\begin{array}{lll}{\sigma }_{\text{ecorr}}^2\hfill & =\hfill & R_X^2({\sigma }_{B_{\text{max}}}^2+{\sigma }_{V_{\text{max}}}^2-2{\sigma }_{[B_{\text{max}},V_{\text{max}}]})\hfill \\ \hfill & \hfill & -2a_X{R}_X({\sigma }_{[\Delta M_{15},V_{\text{max}}]}-{\sigma }_{[\Delta M_{15},B_{\text{max}}]})+R_X^2{\sigma }_{\text{col},\text{intr}}^2.\hfill \end{array}$$(7)

is the uncertainty of the extinction correction and it is asumed that ${\sigma }_{[m_X,B_{\text{max}}]}={\sigma }_{[m_X,V_{\text{max}}]}=0$. Following Marriner et al. (2011) we included a nonzero intrinsic color scatter σ_col,intr, which does not correlate with luminosity. In general, the intrinsic color scatter could be different for each data set, but in our analysis we assumed that it was the same for all data sets σ_col,intr ≃ 0.04 mag. This value was estimated from analysis of the distribution of ${(B-V)}_{B_{\text{max}}}$ as a function of ΔM₁₅ (see Appendix F for details).

The weight factor w in Eq. (5) accounts for the fact that the term σ_ecorr arises only when an extinction correction is applied. Ideally, w is zero for $E{(B-V)}_{\text{obs}}\le 0$ and one for $E{(B-V)}_{\text{obs}}>0$. However, the nonzero intrinsic color scatter makes the situation more complicated. For given stretch factor, among the SNe with $E{(B-V)}_{\text{obs}}<0$ there shouldbe a fraction of intrinsically blue objects that suffered only small reddening. In our analysis these objects will not be extinction corrected, but they should. In an attempt to treat these objects in a more correct way, w was set to linearly decrease from 1 at $E{(B-V)}_{\text{obs}}$ to 0 at $E{(B-V)}_{\text{obs}}=-{\sigma }_{\text{col},\text{intr}}=-0.04$ mag, meaning that an extra uncertainty was assigned to the SNe with small negative $E{(B-V)}_{\text{obs}}$ even though no extinction corrections was applied.

Fig. 7 $V_{B_{\text{max}}}-J_{\text{max}}$ (left) and $V_{B_{\text{max}}}-H_{\text{max}}$ (right) colors vs. $E{(B-V)}_{\text{obs}}$. The dashed, the solid and the dotted lines show the expected relations for Milky Way-type dust and RV = 1.7, 1.9 and 3.1, respectively. Color codes are as in Fig. 6.

4.2.1 NIR Hubble diagram

The model defined by Eqs. (3)–(7) was used to analyze the NIR J and H peak magnitudes of the SNe with z ≥ 0.008. The data was considered consisting of five SN samples: the new SN observations presented in this paper, the BN12 sample, and the published SN sample, which is split into three subsamples - GOLD, SILVER and BRONZE. Two SNe from our sample, SN 2008bz and SNF20080512-008, were observed in the NIR at stretch-corrected phase ~+13 days. This is beyond the range where our fitting method could be used and these two SNe were excluded from the analysis. Further, two SNe from the CSP sample, SNe 2005ku and 2007hx, were found to be significantly brighter than the best-fit Hubble line in H band. CSP (C.Contreras, private communication) confirmed that most of the H band photometry of these two SNe was unreliable. Only one H observation of SN 2005ku was deemed good and the updated photometry was even brighter than the published one. For this reason SNe 2005ku and 2007hx were also excluded. The results of the analysis are shown in Table 11.

The data was first analyzed assuming that SNe Ia are natural standard candles in the NIR and no corrections to the peak magnitudes are needed. The weighted scatter around the Hubble line was 0.154 mag in J and 0.141 mag in H, which is comparable to the best values reported for B band. Inspection of the Hubble residuals revealed a weak dependence on the color excess $E{(B-V)}_{\text{obs}}$ and possibly on ΔM₁₅. When the a_J∕H and R_J∕H correction terms were included in the analysis the scatter of the Hubble residuals only marginally decreased to 0.146 mag in J and 0.136 mag in H. The corresponding H-band Hubble diagram is shown in Fig. 8 with panels (c) and (d) showing the residuals versus Δ M₁₅ and $E{(B-V)}_{\text{obs}}$. The weak dependence of the peak magnitudes on ΔM₁₅ and $E{(B-V)}_{\text{obs}}$, combined with the lack of many highly reddened SNe, is the main reason for the rather small decrease of the scatter. Furthermore, the low values of R_J∕H = 0.42∕0.22 are consistent with Milky Way type extinction with R_V ≃ 1.5 − 2 as commonly found in the optical studies of SNe Ia. The a_X coefficient is ≈ 0.20 in both bands and is smaller than the estimates of Folatelli et al. (2010) and Kattner et al. (2012). The Kattner et al. (2012) analysis shows that the estimate of a_X is highly dependent on the sample used. Thus, the fact that we used a different sample, together with the different approaches to estimate ΔM₁₅ and $E{(B-V)}_{\text{obs}}$, is the most likely the reason for the difference with the previous studies. We also note that the a_X coefficient is much smaller than in the B band (see next section and, e.g., Folatelli et al. 2010), which further supports the conclusion that SNe Ia are better natural standard candles in the NIR than in the optical.

The overall rms around the Hubble line is similar in J and H band, 0.146 ± 0.010 mag and 0.136 ± 0.010 mag, respectively. This is somehow surprising because several previous studies have claimed that H-band Hubble diagram has considerably smaller scatter than J (e.g., Wood-Vasey et al. 2008), but we find that SNe Ia are nearly equally good standard candles in both bands. However, we note that there is no big difference between the J and H band rms’ in the literature samples, but in our and BN12 samples the H band rms is noticeably smaller than in J.

The value of σ_intr is a very important quantity because it tells us how precisely SNe Ia peak magnitudes could be standardized. The uncertainty of the value of σ_ext will affect the estimation of σ_intr from the GOLD, SILVER, and BRONZE samples more compared to BN12 and our samples because the latter have higher mean redshift. The BN12 sample has σ_intr ≃ 0.11 mag in J and even smaller σ_intr ≃ 0.06 mag in H. Even though we assumed rather small σ_ext= 150 km s⁻¹ (300 km s⁻¹ and even 450 km s⁻¹ are sometimes used), the GOLD sample also has σ_intr ≃ 0.1 in both NIR bands. This suggest that the intrinsic luminosity scatter in J and H is most likely around σ_intr ≃ 0.10 mag.

The GOLD, SILVER, and BRONZE have similar redshift distribution and should be equally affected by the uncertainty of σ_ext. Given this, if all sources of uncertainties of the peak magnitudes have been properly accounted for, these three samples should have similar σ_intr. Our analysis shows that σ_intr in the SILVER and BRONZE samples is larger than in the GOLD sample by a factor of about 1.5. Our sample has larger σ_intr in J than any of the other samples, but in H σ_intr is similar to the SILVER sample. This suggest that in spite of our efforts to identify and include all uncertainties that affect the peak magnitudes estimated from incomplete LCs, some additional scatter is still present. One possible source of such scatter is systematic deviation from the s_NIR = s_opt relation. Careful inspection of Fig. 3a reveals that the points in the range s_opt = 0.92–1.02 are systematically above the line, while those in the ranges s_opt = 0.8–0.9 and s_opt = 1.02–1.15 are systematically below it. This raises the interesting possibility that s_NIR has a bi-modal distribution with mean values of s_NIR ≃ 0.8 and s_NIR ≃ 1.05 for the SNe with s_opt below and above ~0.91, respectively, but clearly more observations are needed to confirm whether these systematic trends are real or not. If real, these systematic deviations may affect our and the BRONZE samples differently because the two samples have different stretch distributions. In our sample 4 SNe are in the range s_opt = 0.92–1.04, 11 have s_opt > 1.04 and only one has s_opt < 0.92. In the BRONZE sample on the other hand 8 are in the range s_opt = 0.92 − 1.04, only 5 have s_opt > 1.04 and the majority of 15 SNe have s_opt < 0.92. The SILVER sample has stretch distribution similar to that of the BRONZE and BN12 similar to ours. However, because the SILVER and BN12 SNe have observations close the NIR maximum, these samples should be less sensitive to the uncertainties of the stretch.

Despite the differences of the inferred values of σ_intr and the observed scatter, the SN Ia absolute magnitudes estimated from the GOLD, SILVER, BRONZE and our sample are in a good agreement and differ by lessthan 1σ. This indicates that using under-sampled LCs does not lead to significant bias in estimating the NIR peak magnitudes. However, the BN12 sample is on average 0.15-0.20 mag fainter than the rest. This difference is not negligible and appears to be too large to be explained with differences between the photometric systems used. It is also likely not related to our analysis procedures because Weyant et al. (2014) also noted that the absolute magnitudes inferred from BN12 sample differ from the other samples.

Upon closer examination, the Hubble residuals hinted at possible dependence on the J_max − H_max color (V − [J∕H] were also tested, but the highest degree of correlation was found with the J_max − H_max color). Linear least-squares fit produced nonzero slopes of 0.38 ± 0.11 and − 0.36 ± 0.11 in J and H, respectively.These fits are shown in Fig. 8e and f. The uncertainties were estimated from fitting 10 000 bootstrap replications of the original sample, which with or without preserving the number of SNe in the samples, gave uncertainties of ~0.1. Thus the γ term: $$-{\gamma }_X\times \left[(J_{\text{max}}-H_{\text{max}})-{(J-H)}_0\right],$$(8)

and the relevant error terms were added to Eqs. (3) and (6), respectively. Here ${(J-H)}_0=-0.24$ mag is set to the mean SN Ia J_max − H_max color. The results of the fits with the updates equation also indicate the presence of a significant nonzero γ term. We note that including the γ terms makes the J- and H-band residuals highly correlated and forces γ_J = γ_H + 1. As a results all derived parameters, except for the absolute magnitudes, become the same for J and H (Table 11).

From Table 11, one can see that including the γ term does not lead to decrease of the inferred σ_intr in all samples and in few an increase is observed. However, in all cases the change is within the estimated uncertainties. It can also be seen that in the best observed GOLD and BN12 samples σ_intr does not show a clear decrease and the nonzero γ term is likely not driven by these samples. All this casts doubt of whether the correlation of the Hubble residuals with the J_max − H_max color is real or an artifact of the small sample sizes. However, we note that the range of J_max − H_max color in the GOLD and BN12 samples is smaller than in the others. Thus, the question of whether the Hubble residuals depend on the J_max − H_max color should be answered in the future when a larger sample of well-observed SNe Ia in the NIR becomes available. For this reason we base our further discussion of the results without the γ term.

Fig. 8 H-band Hubble diagram for the full sample with the aH and RH correction terms included. Panels c and d show the Hubble residuals when only aH or RH terms are applied vs. the quantity multiplied by the other, for example, in panel d the Hubble residuals vs. $E{(B-V)}_{\text{obs}}$ when only the aH correction term is applied. Panels e and f show H and J band residuals vs. Jmax − Hmax color. In panel d the dashed lines are the extinction corrections for Milky Way extinction law with RV = 1.7. The two red points surrounded by black circles are SN 2008bz and SNF20080512-008, which have been excluded from the fit (see text). The blue point surrounded a black triangle is PTF10hmv, which was observed only in H.

Fig. 9 B-band Hubble diagram. The right panels show the Hubble residual when only the LC shape or the extinction correction is applied.

4.2.2 B and V bands

The model defined by Eqs. (3)–(7) was also used to analyze the B- and V -band peak magnitudes of the whole nearby SN Ia sample⁹ with redshift above 0.015 plus the 16 SNe from this work – 348 peak magnitudes in total. The term $2R_X({\sigma }_{[B_{\text{max}},V_{\text{max}}]}-{\sigma }_{B_{\text{max}}}^2)$ was added to Eq. (7) because in this case m_X = B_max. The results of the analysis are shown in Table 12 and the B-band Hubble diagram in Fig. 9. First we note that the parameters of our sample fall well in the ranges of the other samples. The overall weighted 1σ scatter around the best-fit Hubble line is 0.176 ± 0.007 mag and 0.172 ± 0.007 mag and B and V bands, respectively.In the individual samples the observed scatter varies between 0.15 and 0.25 mag. The intrinsic luminosity scatter σ_B,intr is between 0.10 and 0.15 mag, and larger by about 0.02 mag in V. Betoule et al. (2014), Scolnic et al. (2014a) and Scolnic et al. (2017) estimate similar B-band intrinsic luminosity scatter for the nearby SN samples, even though using a different LC fitting approaches and applying different cuts to the samples used. From Table 12 one can also see that the absolute magnitudes for the different data sets span a range of ~0.16 mag. Tensions between the nearby SN Ia samples have been reported previously (e.g., Hicken et al. 2009b; Suzuki et al. 2012; Ganeshalingam et al. 2013; Scolnic et al. 2014a) and possibly indicate systematic differences between the data sets, and the peak magnitudes of the SNe in common between the four largest data sets indeed show offsets up to 0.05 mag (Table 5). The nearby SN sample does suffer from selection effects, which can be seen from plots of ${(B-V)}_{B_{\text{max}}}$ and Δ M₁₅ vs. redshift z. Above z ≃ 0.035–0.040 there are no SNe with ${(B-V)}_{B_{\text{max}}}>0.5$ mag, while there are many at the lower redshifts. Similarly, all SNe with z > 0.06 have Δ M₁₅ < 1.4 mag, but at the lower redshifts many SNe with ΔM₁₅ > 1.4 mag are present. However, the distributions of the redshift, ΔM₁₅ and ${(B-V)}_{B_{\text{max}}}$ in the largest samples are very similar and the selection biases are unlikely to be the main cause for the tensions between the sets. Another potential bias is related the differences of the stretch distributions and possibly the luminosity distributions of normal SNe Ia in passive and active host galaxies (see, e.g., Sullivan et al. 2010). The nearby samples contain the same mixture of SNe in passive and active hosts and this bias is probably not important in our case.

Similar to other recent studies we also find low values of R_V = 1.81 ± 0.060 and R_B = 2.90 ± 0.062 even tough intrinsic color scatter of 0.04 mag was included in our analysis. Scolnic et al. (2014b) have claimed that erroneous treatment of the color scatter as luminosity scatter will bias the estimates of a_B and R_B coefficients. R_B in particular can be biased by as much as 30% toward lower values and the authors claim that when the color scatter is treated properly the usual R_B ≃ 4.1 should be recovered. However, from the lower right panel of Fig. 9 one can see that R_B = 4.1, (shown with the dotted line) is inconsistent with the data. Moreover, from analysis of simulated data sets with luminosity and color scatter, and with the same number of SNe and distributions of ${(B-V)}_{B_{\text{max}}}$, Δ M₁₅ and redshift as the full data set, we could only find small bias in R_B at the level of ≤ 10% and even less in a_B. The differences between our results and those of Scolnic et al. (2014b) are likely due to the different color distributions of the SN samples analyzed. The majority of the SNe in the sample simulated by Scolnic et al. (2014b) have ${(B-V)}_{B_{\text{max}}}<+0.2$ mag and none with ${(B-V)}_{B_{\text{max}}}>+0.3$ mag, while the full nearby sample has many SNe with ${(B-V)}_{B_{\text{max}}}=0.2-1.5$ mag. Repeating the analysis only with the SNe with $E{(B-V)}_{\text{obs}}<+0.3$ we find R_B = 3.57 ± 0.13 and R_B = 2.76 ± 0.13 with and without intrinsic color scatter, respectively. Simulated data sets were also analyzed and the results showed that the bias is larger in samples with lower color range. It should also be emphasized that the assumption of universal R_B value for allSNe most likely is not valid for any real SN sample. Indeed, large variations of R_V from event to event have been found by Burns et al. (2014) and Amanullah et al. (2015), with the latter work also including UV measurements, where the differences in extinction laws are the largest.

The slopes of the B_max −ΔM₁₅ and V _max −ΔM₁₅ relations are a_B = 0.737 ± 0.037 and a_V = 0.587 ± 0.037. These values are again lower than those found by Folatelli et al. (2010) and similarly to the NIR, we attribute this to the different sample, and the different approaches we used to estimate Δ M₁₅ and $E{(B-V)}_{\text{obs}}$.

The analysis of the color indexes at maximum light suggested that several SNe may be reddened by dust with R_B ≃ 4.1 and not 2.9(see Fig. 7). The mean color excess of these SNe is ~ +0.3 mag. If these SNe had to be corrected with R_B = 4.1 but we used R_B = 2.9 instead, then their Hubble residuals should be onaverage ~ +0.3 mag. However, this is not the case and the Hubble residuals are on average close to zero. Thus, most likely the behavior of the maximum colors of these SNe is an intrinsic property and not caused by dust extinction with R_B = 4.1. We also find that it is likely to be unrelated to the expansion velocity dichotomy noted by Wang et al. (2009a).

4.3 V − J and V − H color indexes evolution

Krisciunas et al. (2004b) have studied the V −NIR colors of SNe Ia and found that their evolution between phases − 10 and + 10 days is fairly uniform. They also showed that SNe with ΔM₁₅ < 1.1 mag and ΔM₁₅ > 1.1 follow different V − H loci. We studied the evolution of the V − J and V − H colors with the SNe of the GOLD sample with low reddening E(B − V) < 0.2 mag. The observed V, J and H magnitudes, corrected for the Milky Way dust extinction, were converted to rest-frame and the observed epochs were converted to stretch-corrected phases using Eq. (4) with s_opt. The V magnitudes were interpolated at the phases of the NIR observations. Most SNe have well-sampled V band LCs and the interpolation was done with smoothing splines. Only in a few cases interpolation from the templates fits was used.

We assumed that the evolution of the intrinsic colors indexes with the stretch-corrected phase is the same for all SNe. Under this assumption the offsets between color curves of the different SNe arise from different dust extinction in the host galaxies and possibly from intrinsic color scatter. The offset should also be a function of the decline rate parameter ΔM₁₅ because the slope of the “luminosity – LC shape” relation in V and the NIR bands is different. Thus, the observed V − X color curves of each SN were corrected by $${(V-X)}^{\ast }=(V-X)-R_{(V-X)}E(B-V)-a_{(V-X)}(\Delta M_{15}-1.1),$$(9)

where X is either J or H, and R_{(V −X)} = R_V − R_X and a_{(V −X)} = a_V − a_X were computed from the R and a coefficients determined from the Hubble diagram fits. The corrected ${(V-X)}^{\ast }$ curves between phases − 12 and + 12 days were stacked together and weighted least-squares fit with third-order polynomial function was performed. Intrinsic scatter term σ_intr was also added in quadrature to all uncertainties so that the reduced χ² of the fit was one.

After the first iteration it was noted that four SNe had color curves significantly deviating from the rest. SNe 2002dj, 2006le, 2006lf and 2008gp were all too blue at all epochs and were excluded from the analysis. The final sample consisted of 22 LCs of 20 SNe with ΔM₁₅ and $E{(B-V)}_{\text{obs}}$ in the ranges 0.82 −−1.56 mag and − 0.05 to +0.16 mag, respectively. The coefficients of the best-fit polynomial function, the inferred intrinsic scatter and the weighted standard deviation around the fit are shown in Table 13. The corrected color curves up to phases +45 days and the best fits are shown in Fig. 10.

The derived intrinsic scatter of σ_intr ≃ 0.09–0.11 mag is comparable to that reported by Burns et al. (2014), even tough we used one mean value of R_{(V −X)}, while Burns et al. (2014) derived individual values for each SN. We note also that part of the scatter of the V −NIR colors is not statistical. A closer examination of the data shows that most SNe are either bluer or redder than the mean loci at all epochs. This behavior may be related to the use of average R_{(V −X)}, but it may also result from other sources such as systematic differences in photometry and intrinsic color scatter. The fact that four SNe from the low-reddened sample had systematically too blue colors demonstrates this point. We note, however, that the offsets of SNe 2006le and 2006lf are too big to be explained by a different extinction law.

Fig. 10 V − J and V − H un-reddenedloci for a SN with ΔM15 = 1.1. The third-order polynomial fit of the data between phases − 12 and +12 days is also shown. For clarity the error bars are not plotted.

5 Discussion

5.1 SNe Ia in the optical

We analyzed the B and V LCs of our new SNa Ia sample together with the largest nearby SN Ia compilation to date – 458 LCs from 10 different samples. The Hubble diagram of all nearby SNe with redshift above z = 0.015 has weighted scatter around the best-fit line of ~0.176 mag in B and ~0.172 mag in V band. The B band scatterin the individual samples varies (Table 12) and for example it is between 0.15 and 0.25 mag for the four largest samples CSP, Berkeley, CfA3 and CfA4, and our sample has 0.2 mag. Similarly, the corresponding B band intrinsic scatter σ_intr varies in therange 0.1–0.2 mag with our sample having 0.17 mag. These results show that our data fit very well with findings in the literature. The four largest samples have similar redshift distributions and photometric quality, and the differences in σ_intr indicate the likely presence of additional sample-dependent source of uncertainty beyond the intrinsic SN Ia luminosity scatter because the latter should be the same for all samples. Several of the samples that we analized have intrinsic scatter σ_intr = 0.10–0.12 mag. From a similarly heterogeneous low-redshift samples Betoule et al. (2014), Scolnic et al. (2014a) and Scolnic et al. (2017) find an intrinsic luminosity scatter ~0.10–0.12 mag, a value which agrees with our estimate.

Analyzing the distribution of the maximum B − V color of 393 SNe Ia we conclude that the intrinsic B − V color of normal SNe Ia at maximum light is a nonlinear function of Δ M₁₅ and has intrinsic scatter ~0.04 mag. Including this in the analysis, we found that the host galaxy dust extinction is best described with Milky way-type dust with low R_V ≃ 1.8–1.9 in accordance with many other studies. Because of the large span of the colors this estimate is robust and it is easy to see from Fig. 9 that the typical Milky Way value R_V = 3.1 is inconsistent with the data.

Further, restricting the B band Hubble diagam analysis only to the 78 SNe in common with Betoule et al. (2014), which were treated as a single sample, weighted scatter of the Hubble residuals of 0.162 mag and σ_intr = 0.11 mag were obtained. Betoule et al. (2014) selected the best-observed low-reddening normal SNe and our results confirm that such selection criteria provide improvement over the use of all SNe. However, repeating the analysis using only the SNe with $E{(B-V)}_{\text{obs}}<0.3$ mag we find overall scatter of 0.173 mag, thus suggesting that LC quality is a more important contributor to the scatter than the reddening. In general, the similarity of our results to those of Betoule et al. (2014), despite the different analysis approaches, is encouraging given that we only used B and V band LCs.

5.2 SNe Ia in the NIR

The analysis of our new sample of 16 SNe together with the largest compilation of 102 SN Ia NIR LCs to date confirms previous conclusions that SNe Iaare as good, if not better, standard candles in the NIR as in the optical. The overall scatter of the Hubble diagram is ~0.14–0.15 mag in bothJ and H, which is similar to the best published value of 0.14mag in B band (Betoule et al. 2014). The analysis yields R_J ≃ 0.4 and R_H ≃ 0.2, which are consistent with the lowvalue of R_V ≃ 1.8–1.9 derived from the optical observations. This implies that the extinction corrections in the NIR are by a factor of ~7–15 smaller than in the optical B band. The same holds true for the LC shape correction too. We find a_J∕H ≃ 0.20, which is by a factor of 3.5 smaller than in B band. This demonstrates one of the advantages of the NIR observations, namely, that the effect of possible systematic uncertainties associated with the LC shape and the extinction corrections are greatly reduced compared to the optical.

The best-observed nearby GOLD and the BN12 samples indicate rather small intrinsic luminosity scatter σ_intr ≃ 0.10 mag, which is comparable to the latest estimates in the optical B band (Betoule et al. 2014; Scolnic et al. 2014a, 2017). However, one should be cautious about the derived low intrinsic scatter. Estimating this quantity is difficult because it requires identifying all other sources of uncertainty and correctly accounting for them. In addition, low-redshift samples are sensitive to the assumed value of the galaxy velocity dispersion, the exact value of which is rather uncertain. One possible source of uncertainty is the inclusion of the K-correction uncertainties in the error budget. We have assumed that the scatter of the individual K-corrections reflects real difference on the energy distribution in the NIR spectra of different SNe, hence this scatter directly adds to the uncertainty of the peak magnitudes. However, this may not necessarily be the case. It is possible that least part of the observed scatter of the K-corrections is introduced by our procedure to prepare the spectra in the form suitable for computing K-corrections (Appendix E). For example, the linear interpolation across the telluric bands is only a crude approximation. Both the spectra and the photometry to which they were matched may potentially contain calibration errors, and the host galaxy dust extinction estimation may not be accurate enough. All this certainly introduces extra scatter in the individual K-corrections. Adding this extra scatter will lead to deriving smaller σ_intr. Therefore, for accurate determination of the NIR intrinsic luminosity scatter it will be necessary to estimate correctly the uncertainty of the K-corrections. To achieve that, more high-quality photometry and spectroscopy of nearby SNe Ia is clearly needed. Considering this in Tables 6–9 the uncertainties of the peak magnitudes estimated by the LC fitter and the uncertainties of the K-corrections are given separately.

Our findings are also expected from theoretical radiative transfer models. Kasen (2006) found that the NIR peak magnitudes are much less sensitive to the amount of ⁵⁶Ni generated in the explosion (the quantity that to a large extent controls the “luminosity – LC shape” relation) than the optical peak magnitudes. Furthermore, Kasen (2006) found that significant variation of other physical parameters, such as progenitor metallicity, electron capture elements, and production of intermediate mass elements led to only a small scatter of ~0.1 mag of the NIR peak magnitudes, a value which is consistent with the intrinsic scatter found in both optical and NIR.

After applying the LC shape and the extinction corrections it was noticed that the Hubble residuals showed a possible correlation with the J_max − H_max color in both J and H. Adding the J_max − H_max term in the regression analysis leads to further small reduction of the scatter of the whole sample. Even though our analysis indicated that this term is statistically significant, further observations are needed to confirm this finding and whether this correction term will improve the NIR SN Ia standard candle. Recently, the CSP, CfA and SweetSpot programs published updates on their previous data releases and added newly observed objects (Krisciunas et al. 2017; Weyant et al. 2018; Friedman et al. 2015). Most of the new objects are in the Hubble flow and will provide further opportunities to verify the NIR SN Ia standard candle and look for additional correlations such as with the J_max − H_max color.

Recently, the BAO technique has emerged as a very accurate probe for studying the expansion history of the Universe, a very welcome addition to the cosmology toolkit. However, this technique also has certain limitations, most notably the small cosmic volume probed at low redshift (see, e.g., Alam et al. 2017). At low redshift SNe Ia provide significantly better measurements of relative distances and hence finding ways to reduce the systematic uncertainties of the SN Ia technique is particularly relevant. We have seen that rest-frame NIR observations can provide a way to achieve that. Furthermore, SN Ia appear indispensable to constrain models of the dark energy involving time-dependent equation of state parameter (see Figs. 29 and 17 in Anderson et al. (2014) and Alam et al. (2017), respectively). For this, as systematic-free SN Ia observations as possible over a wide redshift range will be needed and in the next section we discuss the prospects for rest-frame NIR observations of SNe Ia at high redshift.

5.3 NIR observations at high redshift

In order to provide constraints of the cosmological parameters with rest-frame NIR observations of SNe Ia, a large sample of events out to redshift z = 1 needs to be obtained. NIR observationsof SNe Ia beyond z ~ 0.5 will be possible with the next generation of 40-m class telescopes (ESO Extremely Large Telescope, The Thirty Meter Telescope) and James Webb Space Telescope (JWST). For example, in the redshift range z = 0.6 − 1.2 the rest-frame J and H bands are shifted to the JWST F200W, F277W and F356W bands, respectively¹⁰. The observed peak (Vega) magnitudes of a typical SN Ia in these bands will be ~23–24 mag and the JWST ETC estimates that with 15 min exposure it should be possible to obtain photometry with S∕N ≃ 20–80 at NIR maximum light and with S∕N ≃ 10–45 at 10 days past $t_{B_{\text{max}}}$. Observations out to z ~ 1.8 will also be possible with exposure times of 40–60 min, but only at maximum light. We have seen that the GOLD sample – the SNe with well-sampled NIR LCs – provides the best estimate of the peak magnitudes and the smallest scatter on the Hubble diagram. However, obtaining LCs with at least four points (plus one more if galaxy reference image is needed) is a considerable investment in terms of observing time at the world’s most advanced and expensive observatory. Thus, the main goal of this study has been to test the accuracy of estimating the NIR peak magnitudes of SNe Ia with only a single observation coupled with good optical data, which is the most efficient approach to use rest-frame NIR observations of SNe Ia in cosmology.

To measure the peak magnitudes from under-sampled LCs, we emplyed a method orgnally used by Krisciunas et al. (2004c,b). We have further verified the applicability of this method by analyzing the optical and NIR LCs of the GOLD sample and discussed in detail the associated uncertainties when under-sampled LCs are used. We have performed series of simulations based on the NIR templates and re-sampling the NIR LCs of the GOLD sample. The results showed that it is possible to estimate the NIR peak magnitudes from a single observation with accuracy better than 0.03 mag if the observation is obtained within a few days from the NIR maximum. If the observation is obtained further away from the maximum the accuracy decreases but is still better than ~0.10 mag.

To test the method on real data we analyzed published LCs together with our new sample of 16 SNe. The published low-redshift sample was divided into three groups GOLD, SILVER, and BRONZE according to the sampling of the NIR LCs, plus the BN12 sample, which is at higher redshift. By design our new sample contains SNe with single NIR observations and only two SNe have two observations. Out of 28 SNe in the BRONZE sample, 16 have at most two observations per band, and all the others have more observations in at least one band. Out of the 15 SNe in the SILVER sample, only four have two or fewer observations. Thus, along with our new sample, the BRONZE one is also well-suited to test the methodology.

The intrinsic luminosity scatter σ_intr of SNe Ia in the NIR was estimated from the GOLD and BN12 samples to be around 0.10 mag in both J and H bands. The scatter inferred from the other samples is larger, most likely because of unaccounted additional uncertainties arising from fitting under-sampled NIR LCs. The intrinsic scatter inferred from our and the BRONZE samples is in the range ~ 0.15 – 0.21 mag and ~ 0.15 – 0.17 mag in J and H, respectively, i.e. larger by a factor of ~2 and ~1.5 than that inferred from the GOLD and BN12 samples. This difference might be related to the much larger flux variation in J than in H after maximum, ~1.3 mag vs. ~0.4 mag, but it may also be an intrinsic property of SNe Ia. Regardless of the reason, our results suggest that when working with under-sampled LCs H band should be preferred over J.

One may then ask the question what is the best observing strategy when limited observing time is available. Should one obtain only well-sampled LCs with four or five points plus a reference image or three times larger sample of SNe with only one observation plus a reference image? Considering the values of the intrinsic luminosity scatter that we estimate in the differentsamples and the results of the simulations, the latter approach may be advantageous, especially if H-band observations can be obtained within two to three days from the NIR maximum.

In the future most SN discoveries at high redshift will likely come from LSST, WFIRST and the DESIRE surveys, which will also have the capability toprovide B and V LCs out to redshift z ~ 1.5. Taking into account that the minimum technical time to trigger JWST will likely be around two days, to observe a SN within a few days from NIR maximum, the SN candidates need to be spectroscopically classified as SNe Ia at least one week before $t_{B_{\text{max}}}$. Very effective selection of SNe Ia for spectroscopic screening (90% success) can be made by studying the early multicolor LCs of the candidates, as shown by Sako et al. (2008). Thus, together with a dedicated program to classify SNe Ia on 8-m class telescopes, LSST and space missions should provide enough targets for NIR follow-up to build NIR Hubble diagram out to redshift z = 1.

6 Conclusions

We obtained multi-epoch optical UBV RI and single-epoch NIR J and H photometry of 16 SNe Ia in the redshift range z = 0.037–0.183. This sample doubles the leverage of the current SN Ia NIR Hubble diagram. We analyzed this new data set together with 102 NIR and 458 optical LCs of normal SNe Ia from the literature. Our findings can be summarized as follows.

A subsample of SNe with well-sampled NIR LCs (GOLD sample) allows template fitting with the stretch and the time of maximum as free parameters. We find that the optical and NIR stretches are nearly one-to-one related and that the time of B-band maximum light estimated from the NIR and optical LC fits are very close. This further reinforces the results of Krisciunas et al. (2004c,b) and Wood-Vasey et al. (2008) who showed that smooth LC templates, the time axis of which is stretched according to the optical stretch parameter, accurately describe the rest-frame J and H-band LCs within 10 days from B-band maximum. This allowed us to estimate the NIR peak magnitudes of the SNe with under-sampled LCs, including our new sample.

Based on simulations, we find that it is possible to estimate the NIR peak magnitudes with accuracy better that 0.025 mag with even a single observation if obtained within two to three days from the time of NIR maximum light. The accuracy decreases to ~0.12 mag if the observation is obtained ten or more days after the time of NIR maximum light.

Our analysis confirms previous findings that SNe Ia are at least as good standard candles in the NIR as in the optical, now with good statistics of >100 objects. We find evidence that the NIR Hubble residuals weakly correlate with the LC shape parameter Δ M₁₅ and the color excess $E{(B-V)}_{\text{obs}}$, and for the first time we report a possible dependence on the J_max − H_max color. The intrinsic luminosity scatter of SNe Ia in the NIR is estimated from the best-observed SNe to be around 0.10 mag. The weighted rms scatter of the Hubble diagram is 0.14 ± 0.01 mag in bothJ and H, and SNe Ia appear to be equally good standard candles in both bands.

The intrinsic luminosity scatter derived from the SNe with under-sampled LCs is larger than that estimated from the best-observed SNe by a factor of ~1.5 in H band and by ~2 in J band. Together with the results of the simulations, this suggests that obtaining three times more SNe with single observation may be advantageous compared to obtaining one SN with well-sampled LC, especially if the observation can be obtained close to the NIR maximum.

Compared to the optical, the corrections for the LC shape in the NIR are by a factor of 3.5 smaller. The extinction corrections are also smaller by an order of magnitude and consistent with low R_V ≃ 1.9. Thus, possible systematic uncertainties associated with the LC shape and the extinction corrections are greatly reduced in the NIR.

We find that the intrinsic B − V color of SNe Ia at B-band maximum light, ${(B-V)}_{B_{\text{max}},0}$, can be described by the following nonlinear function of ΔM₁₅: $${(B-V)}_{B_{\text{max}},0}=0.036\exp \left[\frac{\Delta M_{15}-1.3}{0.3455}\right]-0.084(\pm 0.01).$$

The intrinsic scatter of ${(B-V)}_{B_{\text{max}},0}$ is estimated to be ~0.04 mag.

The weighted scatter of the B band Hubble residuals is 0.176 ± 0.007 mag. The inferred intrinsic luminosity scatter was found to vary between the samples indicating possible systematic uncertainties. The dust extinction in the host galaxies is on average described with R_V ~ 1.8–19. Because of the large span of the maximum B − V color of the nearby SN Ia sample our estimation of R_V is robust and R_V = 3.1 is inconsistent with the data.

In conclusion, we have now extended the NIR SN Ia Hubble diagram to its nonlinear part at z ~ 0.2, that is, where the lever arm in redshift makes it possible to start probing the cosmic constituents, most notably the properties of dark energy. Interestingly, our study confirms that it is feasible to accomplish this result with very modest sampling of the NIR LCs, if complemented by well-sampled optical LCs to determine the SN phase and the optical stretch parameter. This bodes well for the future, in particular when JWST will be operational. JWST will have the capabilities to observe the rest-frame J and H out to redshift z ≃ 1.8, and LSST, WFIRST and DESIRE surveys will provide SNe with well-sampled rest-frame B and V LCs out to redshift z ≃ 1.5. In the view of our results, we suggest that the most efficient way to extend the SN Ia NIR Hubble diagram to high redshifts will be to observe a large samples of SNe Ia with a single observation close to the NIR maximum.

Acknowledgements

This work is based on observations collected at the Nordic Optical Telescope (NOT) located at the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias (La Palma, Spain). Part of the data presented here were obtained in part with ALFOSC, which is provided by the Instituto de Astrofisica de Andalucia (IAA) under a joint agreement with the University of Copenhagen and NOTSA. This paper also uses observations made at the South African Astronomical Observatory (SAAO). Part of the observations were also obtained at the NOT during a student training course in Observational Astronomy provided by Stockholm Observatory and the NorFA Summer School in Observational Astronomy. We are grateful to the NOT astronomers Errki Kankade, Amanda Diupvik, Thomas Augusteijn, Carolin Villforth, Sami Niemi, Tapio Pursimo, and Helena Uthas who performed Service observations, Lison Malo and Celine Reyle for giving up part of their time to observe two of our targets, and Kate Maguire for providing us with optical LCs. We thank the SNfactory team, Greg Aldering, Peter Nugent, and Rollin Thomas in particular, and the SDSS-II team for making their SN discoveries known to us. V.S. acknowledges financial support from Fundação para a Ciência e a Tecnologia (FCT) under program Ciência 2008. This work was partly funded by FCT with the research grant PTDC/CTE-AST/112582/2009. AG and RA acknowledge generous support from the Swedish Research Council. This work has made use of the NASA/IPAC Extragalactic Database (NED), NASA’s Astrophysics Data System, and data products from SDSS and SDSS-II surveys. Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, andthe Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/.

Appendix A Local sequences of stars

The calibrated magnitudes of the local sequences of stars are given in Table A.1 and the stars are indicated on the finding charts in Fig. A.1.

Fig. A.1 Numbers indicate the local sequence stars. In most cases the number is right-up of the star and in a few cases only it is left-down. The SNs are located at the intersection of the two ticks.

Appendix B Second-order atmospheric extinction

The second-order extinction term in U and B bands arisesfrom the fact that the transmission of Earth’s atmosphere is not constant but decreases rapidly at wavelengths below 5000 Å. The change of the transmission across the U and B bands is large enough to cause stars with different spectral energy distributions (SED) to experience slightly different extinction. Thus, the atmospheric extinction in U and B bands is not constant but it is a function of the object SED. It can be parametrized as: $$\begin{array}{lll}{k}_B\hfill & =\hfill & {k^{\prime }}_B+{k^{″}}_B(B-V)\hfill \\ k_U\hfill & =\hfill & {k^{\prime }}_U+{k^{″}}_U(U-B),\hfill \end{array}$$

where ${k^{\prime }}_{U,B}$ and ${k^{″}}_{U,B}$ are the first-order and the second-order extinction. The second-order extinction can be neglected in V RI. Unfortunately, in the majority of the photometric works on SNe Ia the second-order extinction term has been neglected. Here we discuss the uncertainties that are caused by this.

The usual practice in SN photometry is to calibrate the magnitudes of the stars in the field of the SN in several photometric nights. This is usually done by deriving the parameters in Eqs. (1) by observing Landolt (1992, 2009) standard fields that contain both blue and red stars in photometric conditions, and applying the equations to the field stars. Then in nonphotometric nights the SN magnitudes are calibrated using the field stars as tertiary standards. Simple rearrangement of Eqs. (1), taking B band as an example gives: $$B=b+(c{t}_B-{k^{″}}_{B}X)(B-V)-{k^{\prime }}_{B}X+z{p}_B,$$(B.3)

If ${k^{″}}_B$ is not fitted (or fixed to a reasonable value ~−0.03), one would not measure the true color term ct_B but rather $c{t}_B-{k^{″}}_{B}X$, meaning that the measured color-term will depend on the airmass of the observations. The second-order extinction in U and B has been measured and it is about − 0.03. Using the spectrophotometric atlas of Stritzinger et al. (2005), various U and B filters, and atmospheric transmission measured at different observatories we estimated ${k^{″}}_B$ and ${k^{″}}_U$ by means of synthetic photometry to be ~−0.03 as well. Therefore, two standard fields observed airmass X = 1 and 2 will give color terms that differ by ~0.03.

Another problem arises because the B − V color index of SNe Ia changes from ~−0.1 around maximum light to 1.3 about a month later. If a SN is observed at airmass 1.5 the correction for the second-order extinction at maximum and at +30 days will be ~0.0 and 0.058 mag, respectively. When the second-order extinction term is not included in the calibration, this will not be corrected for. While around maximum the effect will be negligible, one to two months after maximum the B magnitudes maybe as much as 0.06 mag in error. This will directly affect the B − V color index, which is often used to estimate the host galaxy extinction. Therefore, the estimation of the peak magnitudes can be affected by as much as 0.18 mag for Milky Way type dust with R_V = 3.1.

Appendix C S-corrections

C.1 Photometric systems responses

To recover the effective system responses we follow the methods described in Stanishev et al. (2007). The relevant information was either taken from Stanishev et al. (2007) or from the corresponding observatories. Whenever possible measured filter transmission and CCD QEs were used; if unavailable, specifications from the manufacturer were used. The total system responses were computed by multiplying the filter transmissions by (a) the CCD quantum efficiency (QE), (b) the reflectivity of at least two aluminum surfaces, (a) a telluric absorption spectrum at airmass 1.3. The lens and window transmissions were not included because this information was unavailable. We note that the transmission of the Earth atmosphere was not included in the calculation of the effective U and B bands. The reason is that when computing the instrument color terms the second-order extinction has been taken into account. The effect of the atmospheric transmission on the U and B bands is to shift the bands effective wavelength redward with increasing the airmass of the observations. Taking into accountthe second-order extinction meant that the computed color terms correspond to observations above the atmosphere.

Synthetic magnitudes were calculated from Stritzinger et al. (2005) spectrophotometry of Landolt standard stars, $$m^{\text{syn}}=-2.5\log \left({\displaystyle \int }{f}_{\lambda }^{\text{phot}}(\lambda )R^{\text{nat}}(\lambda )d\lambda \right).$$(C.1)

The difference between the synthetic and the observed photometry was fitted as a function of the observed color indices to compute synthetic color-terms (ct^syn), for example, for B we have $$B^{\text{std}}-B^{\text{syn}}=c{t}^{\text{syn}}(B^{\text{std}}-V^{\text{std}})+z{p}^{R^{\text{nat}}}.$$(C.2)

For the V RI-bands, the ct^syn’s were close to the observed ones ct^obs. In some cases small differences exceeding the uncertainty were accounted for by shifting the filter transmissions until ct^syn matched ct^osb. $z{p}^{R^{\text{nat}}}$ is the ZP for synthetic photometry with the band R^nat.

As in Stanishev et al. (2007) we found that the synthetic U and B bands were always too blue. This is, to some extent, to be expected because the neglected optical elements such as lenses or windows, anti-reflection and other coatings will tend to reduce the system sensitivity shortward of ~4000 Å. To account for the net effect of these uncertainties we modified the U and B bands by multiplying them with a Sigmoid function: $$F(\lambda ;{\lambda }_0,\Delta )=\frac{1}{1+\exp (-(\lambda -{\lambda }_0)/\Delta )}.$$(C.3)

The parameters λ₀ and Δ were selected so that ct^syn matched ct^obs for both U and B.

Following Stanishev et al. (2007) we also rederived the standard Johnson-Cousins system responses with the revised filter functions of Bessell & Murphy (2012). This procedure also gave the ZPs $z{p}^{R^{\text{std}}}$ for synthetic photometry with the standard system responses R^std.

C.2 Computing the S-corrections

The procedure that we used to calculate the S-corrections was the same as that in Stanishev et al. (2007), but slightly modified to account for the fact that the U and B bands did not include the atmospheric transmission. The calibration correction from natural magnitudes at airmass X, $m_X^{\text{nat}}$, to standard magnitudes m^std is computedby means of synthetic photometry $$m^{\text{std}}=m^{\text{nat}}+S,$$(C.4)

where S is the S-correction: $$\begin{array}{lll}S\hfill & =\hfill & -2.5\log \left({\displaystyle \int }{f}_{\lambda }^{\text{phot}}(\lambda )R^{\text{std}}(\lambda )d\lambda \right)+z{p}^{R^{\text{std}}}\hfill \\ \hfill & \hfill & +2.5\log \left({\displaystyle \int }{f}_{\lambda }^{\text{phot}}(\lambda )R^{\text{nat}}(\lambda )p{(\lambda )}^{X}d\lambda \right)-z{p}_X^{R^{\text{nat}}}.\hfill \end{array}$$(C.5)

Here $f_{\lambda }^{\text{phot}}(\lambda )$ is the photon spectral energy distributions (SED) of the SN, R^nat(λ) and R^std (λ) are the responses of the natural and the standard systems, respectively, and p(λ) is the transmission of the atmosphere at airmass X = 1 above the site of observation. Here $z{p}_X^{R^{\text{nat}}}$ is the ZP obtained by fitting Eq. (C.2) to the synthetic magnitudes $m_X^{\text{syn}}$ computed with the natural passband multiplied by the Earth’s atmosphere transmission at airmass X, for example $$m_X^{\text{syn}}=-2.5\log \left({\displaystyle \int }{f}_{\lambda }^{\text{phot}}(\lambda )R^{\text{nat}}(\lambda )p{(\lambda )}^{X}d\lambda \right).$$(C.6)

Comparing Eq. (B.3) and Eq. (C.4) we see that the term $c{t}_B(B-V)-{k^{″}}_{B}X(B-V)$ is simply an approximation of the S-correction: $$S\simeq c{t}_B(B-V)-{k^{″}}_{B}X(B-V).$$(C.7)

By the way the ZPs in Eq. (C.5) are derived, S ≃ 0.0 for A0 V stars with all color indexes zero, for example, Vega. This ensures that for normal stars the synthetic S-correction will work the same way as the linear color-term corrections.

The best way to compute the S-corrections for a SN is to use spectrophotometry of the same object. If that is not possible the mean spectrum of SNe Ia should be used. In our work we used Hsiao et al. (2007) spectral template, corrected for the Milky Way reddening and the redshift of each SN.

Appendix D s_opt −ΔM₁₅ relation

To derive the new s_opt −ΔM₁₅ relation for normal SNe Ia a subsample of ~200 SNe with well-sampled LCs and observations within two days from $t_{B_{\text{max}}}$ was used. The observed B-band LC were corrected to rest-frame B and normalized to zero peak magnitude. The LCs were divided into 0.05 stretch bins and let s_c be the stretch in the middle of a given bin. The time axes of the LCs that fall into given bin were converted to phases, corrected to the stretch in the middle of the bin by $$\varphi =s_{\text{c}}\frac{t-t_{B_{\text{max}}}}{s_{\text{opt}}(1+z)},$$(D.1)

and stacked together. Smoothing spline was then used to estimate the time maximum and the Δ M₁₅ parameter. The resulting s_opt −ΔM₁₅ pairs are shown in Fig. D.1 together with the third-order polynomial fit (the solid line). It is seen that our relation is close to the one derived by simply stretching the B-band LC template, indicating that the stretched template describes the LC of normal SNe Ia rather well. On the other hand the linear relation derived by Jha et al. (2006), which has been used by many authors, is applicable only in the limited stretch range 0.85–1.05. The parameters of the third-order polynomials to convert from s_opt to Δ M₁₅ and vice versa are given in Table D.1.

Fig. D.1 sopt −ΔM15 relation from this work.

Appendix E NIR K-corrections

To compute the NIR K-corrections we used spectrophotometry of nearby SNe Ia that also had NIR photometry available (Table E.1). For the majority of the available spectra the J,H and K parts were obtained separately with different, nonoverlapping instrument settings. While the relative flux calibration of each individual part is probably accurate enough, the cross-band flux calibration may not be. Besides, the flux across the strong telluric absorption bands at ~1.38 μm and ~1.88 μm was not available. Therefore, before computing the K-corrections the spectra had to be “prepared”. The following procedure was applied.

The noisy parts at the edges of the telluric bands were removed and the fluxacross the bands was linearly interpolated. Next, JHK synthetic magnitudes were computed with the appropriate effective system responses¹¹ and were compared with the observed ones. When necessary the J, H, and K parts were scaled, the interpolation across the telluric bands was repeated and new synthetic photometry was computed¹². The procedure was iterated until the synthetic and the observed colors matched. In the cases when the whole NIR spectral range was observed simultaneously we always found a good match between the synthetic and the observed colors, and only interpolation across the telluric bands was necessary. The spectra were corrected for dust reddening in the Milky Way and the SN host galaxies using E(B − V) values from the Galactic dust map (Schlegel et al. 1998) and from the corresponding photometric papers, and corrected to rest-frame wavelengths. The spectra at similar phases were plotted together and examined. Several spectra with poor S/N, inadequate interpolation across the telluric bands or other problems most likely related to the flux calibration were identified and discarded. Details of the final set of spectra are given in Table E.1. We note that a few spectra do not cover all JHK bands. In addition, a few spectra had issues in one band, but were of sufficient quality in the others.

Fig. E.1 Example of computing the NIR J and H K-corrections at z = 0.02 (left panels) and z = 0.07 (right panels).

Fig. E.2 J and H filters responses used to calculate the NIR K-corrections.

To compute K-corrections at any SN phase smoothing splines were used to interpolate the K-corrections computed from the individual spectra. The smoothing splines havegreat flexibility to describe the data to any pre-defined rms around the fit. We have computed J − J and H − H K-corrections for several redshifts between z = 0.002 and z = 0.2 and visuallyexamined the interpolation to tabulate the optimal smoothing parameter as a function of the redshift. This is a somehow subjective procedure and we tried to select thesmoothing parameters in such a way that the data was well-fitted, but with special care to avoid over-fitting. Examples of the outcome of this procedure are shown in Fig. E.1 for two redshifts of z = 0.02 and z = 0.07. The scatter around the smooth curves increases with the redshift and also depends on partucular filters use (see below). The scatter at redshift z = 0.02 is ~0.01 mag in all bands, at z = 0.08 it increases to ~0.033 and 0.04 mag, for J and H, respectively, and at z = 0.14 the J and H scatter is ~ 0.05 and 0.08 mag. These values and the observed increase with the redshift are similar to the results of Boldt et al. (2014). We note however, that we did not color-correct the spectra to follow single color curve as in Hsiao et al. (2007) and Boldt et al. (2014). As a comparison in Fig. E.1 are also shown the K-corrections computed from Hsiao et al. (2007) template. One can see that there are considerable differences for some SN phases, including before +10 days, which is of interest for our study. Boldt et al. (2014) argue that the dominating part of the scatter is due to intrinsic differences in the energy distribution of different SNe. Thus, the scatter of the individual K-corrections around the smoothing spline function between phases −10 and +10 days was tabulated as a function of redshift and added in quadrature to the uncertainties of the peak magnitudes.

Fig. E.3 Panel a: ${(B-V)}_{B_{\text{max}}}$ color of the SNe used in this work as function of Δ M15. The solid line shows the relation defined by Eq. (F.2). The dashed line has a slope of 0.12 and is shown for a comparison; panel b: Histograms of ${(B-V)}_{B_{\text{max}}}-{(B-V)}_{B_{\text{max}},0}$ for SNe with ΔM15 < 1.1 (182 objects), 1.1 < ΔM15 < 1.3 (163 objects) and ΔM15 > 1.5 (48 objects) with the corresponding fits with Eq. (F.1).

There is no standard photometric system in the NIR, and often slightly different filter sets are employed by different observatories. In our analysis we decided to apply K-corrections from the different observed filter systems to a single JHKs photometricsystem, which we chose to be the one defined by the Mauna Kea Observatories Near-Infrared Filter Set (MKO; Tokunaga et al. 2002). As a practical realization of the MKOsystem we selected the filter and detector combinationused by the NOTCam instrument at the NOT. The observed filter systems that were used in this work are

• (i)

  the MKO filters for the 16 new SNe from this work and the SNe from BN12.

• (ii)

  2MASS (Cohen et al. 2003) for the CfA data set (Wood-Vasey et al. 2008)

• (iii)

  the filters used by The Carnegie Supernova Project (CSP), which are similar to Persson et al. (1998). The CSP filters were used for the CSP data set (Contreras et al. 2010; Stritzinger et al. 2011), the data set from K. Krisciunas and collaborators (Krisciunas et al. 2004c,b, 2003, 2001, 2000, 2007) (the Kri set) and the following SNe from various sources (the Var set) SN 1998bu (Hernandez et al. 2000; Mayya et al. 1998; Jha et al. 1999), SN 1999ee (Hamuy et al. 2002), SN 2000E (Valentini et al. 2003), SN 2002dj (Pignata et al. 2008), SN 2003cg (Elias-Rosa et al. 2006), SN 2003du (Stanishev et al. 2007), SN 2003hv (Leloudas et al. 2009), SN 2004eo (Pastorello et al. 2007a), SN 2005cf (Wang et al. 2009b; Pastorello et al. 2007b), and SN 2008Q (Stanishev et al., in preparation).

• (iv)

  the filters of WHIRC instrument were used for SN 2011fe (Matheson et al. 2012).

and the employed MKO, 2MASS and CSP filter responses are shown in Fig. E.2. K-corrections from MKO to MKO, CSP and 2MASS for several redshifts and phases between − 16 and + 40 days are available at the CDS.

Strictly speaking the above procedure will work accurately only if the observations are obtained with filters closely resembling those of the photometric system to which the photometry is tied. This is the case only for the CSP, CfA and BN12 data sets, and to a large extent for the Kri data set. The observations of the other SNe were obtained with filters that differ from those of the photometric systems to which the photometry was calibrated. Most of our new photometry is tied to 2MASS system. The rest of the SN photometry is tied to Persson et al. (1998) or Hunt et al. (1998) systems, which are similar to each other. Because of this mismatch, it is expected that the photometry will contain small systematic uncertainties. To estimate the magnitude of this uncertainty we computed synthetic photometry with the filters of the four photometric systems considered in this work at two epochs from the Hsiao spectral template – at NIR maxima and a week after $t_{B_{\text{max}}}$. We found that the errors should be smaller than 0.02 mag. The only exception was H band a week after $t_{B_{\text{max}}}$. In this case the error could be as big as 0.06 mag for most SNe and only for WHIRC observations of SN 2011fe it could reach 0.1 mag.

Appendix F (B − V) color at B band maximum

Phillips et al. (1999) were the first to suggest that the intrinsic B − V color of SNe Ia at maximum ${(B-V)}_{B_{\text{max},0}}$ is a linear function of ΔM₁₅ with a slope of ~0.11. The parameters of the ${(B-V)}_{B_{\text{max}},0}-\Delta M_{15}$ relation were determined by linear fit of the bluest objects (see also Folatelli et al. 2010). The disadvantage of this approach is that it requires a prior knowledge of which SNe are likely to be unreddened and that only a small fraction of all available data is used. Here we followed a different path to study the distribution of observed ${(B-V)}_{B_{\text{max}}}$ colors with the much larger sample available today. Nobili & Goobar (2008) have shown that the intrinsic SN Ia ${(B-V)}_{B_{\text{max}}}$ color has a dispersion of about 0.05 mag. The observed colors may be further reddened by interstellar dust along the line of sight, evidence for which is clearly seen as an extended red tail of the color distribution in the nearby SN Ia sample. Let us assume as in Jha et al. (2007) that (i) the intrinsic ${(B-V)}_{B_{\text{max}}}$ colors follow normal distribution with mean c₀ and standard deviation σ_c, and (ii) that the reddening E(B − V) follows exponential distribution with scale τ_e for E(B − V) ≥ 0 mag and is zero otherwise. Under these assumptions the probability density function (PDF) of the observed ${(B-V)}_{B_{\text{max}}}$ color, $\rho ({(B-V)}_{B_{\text{max}}};c_0,{\sigma }_{\text{c}},{\tau }_{\text{e}})$, is the convolution of the two functions (Jha et al. 2007; Marriner et al. 2011): $$\rho (c;c_0,{\sigma }_{\text{c}},{\tau }_{\text{e}})\propto {\displaystyle {\int }_0^{+\infty }\exp }\left(-\frac{{(c-c_0-\tau )}^2}{2{\sigma }_{\text{c}}^2}\right)\exp \left(-\frac{\tau }{{\tau }_{\text{e}}}\right)d\tau ,$$(F.1)

where c is the observed ${(B-V)}_{B_{\text{max}}}$ color and the integration is over the color excess from extinction, τ.

Figure E.3a shows the ${(B-V)}_{B_{\text{max}}}$ colors measured with uncertainty less than 0.04 mag from 393 LCs. The dotted line has a slope of 0.12. It describes the blue edge of the distribution of the SNe with Δ M₁₅ ≤ 1.4 well, but beyond that the colors are clearly increasingly redder than the line. It is clear that a nonlinear function is need to describe the average intrinsic ${(B-V)}_{B_{\text{max}}}$ color as a function of ΔM₁₅. We tried different functions and found that a simple exponential function worked best, while a second-order polynomial was not flexible enough to describe the data.

The parameters of the exponential function that describes ${(B-V)}_{B_{\text{max}},0}$ were determined by fitting Eq. (F.1) to the histogram of $E{(B-V)}_{\text{obs}}={(B-V)}_{B_{\text{max}}}-{(B-V)}_{B_{\text{max}},0}$ computed in 0.05 mag wide bins. Unfortunately, there is no unique way to define the quantity to be minimized. We tried several approaches, which are appropriate under the reasonable assumption that the intrinsic scatter σ_c is independent on ΔM₁₅. First the histogram of the whole data set was fitted and the parameters that give the lowest σ_c was selected. The reasoning behind this is that with incorrect parameters the histogram will be additionally blurred and hence a larger σ_c will be derived; the correct parameters should thus give the lowest σ_c. However, this approach was unsuccessful and the blue edge of the distribution of the SNe with Δ M₁₅ ≥ 1.4 was not described well. The likely reason for this is that the number of SNe with ΔM₁₅ > 1.4 is small and the final solution is mostly driven by the SNe with Δ M₁₅ < 1.4. The approach that we found to work best to address this problem was to split the sample into three bins – SNe with Δ M₁₅ < 1.1 (182 points), 1.1 ≤ΔM₁₅ < 1.5 (163 points) and ΔM₁₅ ≥ 1.5 (48 points) – and fit the histogram of $E{(B-V)}_{\text{obs}}$ in each bin independently. In the fitting the c₀ parameter in Eq. (F.1) was treated as an unknown constant. It is easy to see that with the correct ${(B-V)}_{B_{\text{max}},0}-\Delta M_{15}$ relation, σ_c should be the same in all three bins and the c₀’s should be close to zero. Thus, the quantity to minimize was the scatter of σ_c around the mean in the three bins and the scatter of c₀ around zero added in quadrature. The τ_e parameters were not required to be the same in the three bins because the distribution of the reddening may be different. For example, the SNe with larger ΔM₁₅ are fainter and this may lead to bias against discovering reddened SNe.

With the above procedure, the following function was derived $${(B-V)}_{B_{\text{max}},0}=0.036\exp \left[\frac{\Delta M_{15}-1.3}{0.3455}\right]-0.084(\pm 0.01).$$(F.2)

This relation is shown with the solid line in Fig. E.3a and the histograms in the three bins and the best fits in Fig. E.3b. The mean value of σ_c in the three bins is 0.045 mag with estimated uncertainty of 0.01 mag. The uncertainty of the free term in Eq. (F.2), ± 0.01 mag, was estimated from the fits of the histograms, which give uncertainty of ~0.01 mag of c₀. The constant 1.3 was derived empirically. Several values between 1.1 and 1.4 were tried and the one providing the solution with smallest σ_c was selected. Regarding τ_e, the following values were derived in the three bins, τ_e = 0.14 ± 0.02 mag, τ_e = 0.13 ± 0.02 mag and τ_e = 0.11 ± 0.04 mag. We note that the results are somehow dependent on the binning on ΔM₁₅ and the color, and color error cut (0.04 mag in our case). When a larger sample of SNe,

especially with ΔM₁₅ > 1.4, becomes available it will be possible to estimate the parameters more accurately.

The values of σ_c and τ_e that we findare close to those of Jha et al. (2007), Ganeshalingam et al. (2010) and Folatelli et al. (2010), which were estimated from the late-time SN Ia colors at + 35 days using the Lira-relation (Lira 1995; Phillips et al. 1999). The Lira-relation is considered as one of the best methods to estimate the host galaxy reddening of SNe Ia. In another analysis of late-time B − V color curves of SNe Ia Burns et al. (2014) introduced a new stretch parameter of the B − V curves, s_BV. The authors estimated the intrinsic scatter of the B − V color to be 0.06−0.09 mag. The intrinsic scatter of ${(B-V)}_{B_{\text{max}}}$ that we estimate is smaller than the scatter of B − V at late times derived by the above authors. This indicates that at least for normal SNe Ia the ${(B-V)}_{B_{\text{max}}}$ color shouldprovide as good estimate of the reddening as the Lira-relation. This may be important for high-redshift studies where the observations do not extend to late phases and the extinction must be estimated from the colors at maximum.

As a final note, we would like to point out that the width of the observed color distribution is not determined only by the intrinsic color scatter σ_c and the distribution of the reddening (here parametrized by τ_e). Several other other factors, such as the observational errors, the fact that S-corrections were not applied to all data and possible deviation of the fitting LC templates from the real SN LCs will all increase the width of the observed distribution and lead to overestimate of σ_c and τ_e. Some of these sources of errors may be different for the different samples and therefore σ_c should be regarded as a mean value for the whole nearby SN Ia sample. The mean uncertainty of ${(B-V)}_{B_{\text{max}}}$ in the sample used in this analysis is ~ 0.015 mag and taking it into account the true intrinsic scatter σ_c is likely ~ 0.04 mag.

References

All Tables

All Figures

Fig. 1 griJH and UBV RIJH LCs and the best fits. For the three SDSS-II SNe (2007hz, 2007ih, 2007ik) from bottom to top are plotted g, r − 0.4, i − 1, J − 2.8 and H − 3.8 (no J for SN 2007ik). For the other 13 SNe are plotted U + 0.8, B, V − 0.2, R − 0.6, I − 1.6, J − 3 and H − 3.5. We note that the simultaneous gri or UBV RI fits were only used for these plots. In the analysis we use the results from the simultaneous fits of the g and r, and B and V bands only.

In the text

	Fig. 2 Data used to create the J and H LC templates. The templates areoverplotted with solid lines. σw is the weighted rms of the data around the template during the phases earlier than +10 days.
In the text

	Fig. 3 Panel a: NIR vs. optical stretch parameter. The solid line is the best-line fit. The four encircled points are the SNe shown in Fig. 4; Panel b: histogram of the difference between time of B band maximum estimated from the NIR and optical LC fits from the GOLD sample.
In the text

	Fig. 4 Example of NIR LC fit of SNe from the GOLD sample. The solid line is the best fit, and dotted lines are fits with sNIR = sopt. The J band is shown in blue and H in red.
In the text

	Fig. 5 Peak magnitude uncertainties of a SN with sNIR = 1.0 vs. the phase of the single observation. The 1σ uncertainties of the stretch and the time of maximum were assumed to be 0.05 and 0.44 days, respectively.
In the text

	Fig. 6 Panel a: Jmax − Hmax color vs. ΔM15. The solid line shows the least-squares linear fit to the GOLD sample; panel b: Jmax − Hmax $vs.E{(B-V)}_{\text{obs}}$. The solid line is the least-squares linear fit to the GOLD sample, the dashed and the dotted lines are the expected relation for Milky Way-type dust and RV = 1.7 and 3.1, respectively.
In the text

	Fig. 7 $V_{B_{\text{max}}}-J_{\text{max}}$ (left) and $V_{B_{\text{max}}}-H_{\text{max}}$ (right) colors vs. $E{(B-V)}_{\text{obs}}$. The dashed, the solid and the dotted lines show the expected relations for Milky Way-type dust and RV = 1.7, 1.9 and 3.1, respectively. Color codes are as in Fig. 6.
In the text

Fig. 8 H-band Hubble diagram for the full sample with the aH and RH correction terms included. Panels c and d show the Hubble residuals when only aH or RH terms are applied vs. the quantity multiplied by the other, for example, in panel d the Hubble residuals vs. $E{(B-V)}_{\text{obs}}$ when only the aH correction term is applied. Panels e and f show H and J band residuals vs. Jmax − Hmax color. In panel d the dashed lines are the extinction corrections for Milky Way extinction law with RV = 1.7. The two red points surrounded by black circles are SN 2008bz and SNF20080512-008, which have been excluded from the fit (see text). The blue point surrounded a black triangle is PTF10hmv, which was observed only in H.

In the text

	Fig. 9 B-band Hubble diagram. The right panels show the Hubble residual when only the LC shape or the extinction correction is applied.
In the text

	Fig. 10 V − J and V − H un-reddenedloci for a SN with ΔM15 = 1.1. The third-order polynomial fit of the data between phases − 12 and +12 days is also shown. For clarity the error bars are not plotted.
In the text

	Fig. A.1 Numbers indicate the local sequence stars. In most cases the number is right-up of the star and in a few cases only it is left-down. The SNs are located at the intersection of the two ticks.
In the text

	Fig. D.1 sopt −ΔM15 relation from this work.
In the text

	Fig. E.1 Example of computing the NIR J and H K-corrections at z = 0.02 (left panels) and z = 0.07 (right panels).
In the text

	Fig. E.2 J and H filters responses used to calculate the NIR K-corrections.
In the text

Fig. E.3 Panel a: ${(B-V)}_{B_{\text{max}}}$ color of the SNe used in this work as function of Δ M15. The solid line shows the relation defined by Eq. (F.2). The dashed line has a slope of 0.12 and is shown for a comparison; panel b: Histograms of ${(B-V)}_{B_{\text{max}}}-{(B-V)}_{B_{\text{max}},0}$ for SNe with ΔM15 < 1.1 (182 objects), 1.1 < ΔM15 < 1.3 (163 objects) and ΔM15 > 1.5 (48 objects) with the corresponding fits with Eq. (F.1).

In the text