Type II Supernovae as Distance Indicators at Near-IR Wavelengths

Motivated by the advantages of observing at near-IR wavelengths, we investigate Type II supernovae (SNe II) as distance indicators at those wavelengths through the Photospheric Magnitude Method (PMM). For the analysis, we use $BVIJH$ photometry and optical spectroscopy of 24 SNe II during the photospheric phase. To correct photometry for extinction and redshift effects, we compute total-to-selective broadband extinction ratios and $K$-corrections up to $z=0.032$. To estimate host galaxy colour excesses, we use the colour-colour curve method with the $V\!-\!I$ versus $B\!-\!V$ as colour combination. We calibrate the PMM using four SNe II in galaxies having Tip of the Red Giant Branch distances. Among our 24 SNe II, nine are at $cz>2000$ km s$^{-1}$, which we use to construct Hubble diagrams (HDs). To further explore the PMM distance precision, we include into HDs the four SNe used for calibration and other two in galaxies with Cepheid and SN Ia distances. With a set of 15 SNe II we obtain a HD rms of 0.13 mag for the $J$-band, which compares to the rms of 0.15-0.26 mag for optical bands. This reflects the benefits of measuring PMM distances with near-IR instead of optical photometry. With the evidence we have, we can set the PMM distance precision with $J$-band below 10 per cent with a confidence level of 99 per cent.


INTRODUCTION
Type II supernovae (SNe II) are the explosive end of massive stars (MZAMS > 8 M ⊙ ) that retain an important amount of hydrogen in their envelopes at the moment of the explosion. These events, consequence of the gravitational collapse of their iron cores, are characterized by a luminosity comparable to the total luminosity of their host galaxies, which make them interesting objects for distance measurements. ⋆ This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile. † E-mail: olrodrig@gmail.com The pioneering work of Kirshner & Kwan (1974) marks the beginning of the use of SNe II as distance indicators. In their work they applied the Expanding Photosphere Method (EPM, a variant of the Baade-Wesselink method) to two SNe II, using optical photometry and spectroscopy during the photospheric phase (the phase between the maximum light and the transition to the radioactive tail) to estimate angular and physical sizes, respectively. For the first implementation of the EPM, Kirshner & Kwan (1974) assumed SNe II emit like blackbodies. Years after, Wagoner (1981) demonstrated that the flux of SNe II is diluted as a consequence of their scattering-dominated atmospheres, making necessary SN II atmosphere models to quantify that effect and thus c 2018 The Authors to correct derived distances. Since then, the EPM has been applied using different theoretical atmosphere models (e.g., Eastman et al. 1996;Dessart & Hillier 2005) to an ample number of SNe II (e.g., Schmidt et al. 1992Schmidt et al. , 1994Dessart et al. 2008;Jones et al. 2009;Bose & Kumar 2014;Gall et al. 2016Gall et al. , 2018, where the typical EPM distance precision is found to be about 15 per cent. Empirically, Hamuy & Pinto (2002) found a correlation between the bolometric luminosity at 50 d since explosion and the expansion velocity of the photosphere at the same epoch. This is due to the fact that a more energetic explosion corresponds to a more luminous SN with higher envelope expansion velocities. The latter correlation is the basis of the Standardized Candle Method (SCM), which allows to estimate distances using photometry and expansion velocities inter-or extrapolated at 50 d since explosion. The SCM has been applied to several SN II sets (e.g., Nugent et al. 2006;Poznanski et al. 2009;Olivares E. et al. 2010;D'Andrea et al. 2010;de Jaeger et al. 2015de Jaeger et al. , 2017bGall et al. 2018), yielding a distance precision about 12-14 per cent.
Despite apparent differences between the EPM and the SCM, Kasen & Woosley (2009) showed that the SCM is a recasting of the EPM at 50 d since explosion. Additionally, by means of SN II models, they proposed a generalization of the SCM, which can be applied in any epoch during the photospheric phase. The same idea was investigated empirically by Rodríguez et al. (2014, hereafter R14), who called it the Photospheric Magnitude Method (PMM) to the SCM generalization. Measuring distances with all expansion velocities available during the photospheric phase, decreases observational errors and reduces uncertainties introduced by the interpolation/extrapolation at a certain fiducial epoch. The PMM distance precision is around 6-11 per cent (R14).
For the EPM, SCM, and PMM, optical spectroscopy is necessary in order to estimate expansion velocities. Since the spectroscopy is more time consuming than photometry, expansion velocities are not always available. For this reason, de Jaeger et al. (2015) proposed a method based solely on photometry to standardize SNe II, known as the Photometric Colour Method (PCM). de Jaeger et al. (2017b) applied the PCM to a SN II sample with redshift up to 0.5, finding that the PCM distance precision is around 17 per cent.
Most of distances measurements with the latter methods have been performed with optical photometry. However, observing at near-IR wavelengths has two clear benefits that in principle can improve their use as distance indicators: 1. Near-IR light is less affected by dust. Methods to measure colour excess due to SNe II host galaxies (e.g., Schmidt et al. 1992;Krisciunas et al. 2009;Olivares E. et al. 2010;Poznanski et al. 2012;R14;Pejcha & Prieto 2015) are still not well established. Therefore, it is propitious to observe at near-IR wavelengths in order to reduce the effect of miscalculation of the colour excess. Moreover, the estimation of a representative extinction curve along the SN II line of sight is still controversial. Assuming the family of extinction curves of Cardelli et al. (1989), some studies are in favour of a Galactic RV = 3.1 (e.g., Pejcha & Prieto 2015), while other authors found results in favour of lower values (Poznanski et al. 2009;Olivares E. et al. 2010;de Jaeger et al. 2015). Since the choice of a certain extinction curve has more impact at optical than at near-IR wavelengths (e.g., Schlafly et al. 2016), it is preferable to perform photometric observations at those wavelengths in order to diminish systematics induced by the assumption of an incorrect extinction curve.
2. Contamination by metal lines is less severe at near-IR wavelengths. Among the few metal lines identified in the near-IR, we remind: in the J-band range there is a feature at λ = 1.2 µm possibly due to a Si i multiplet (Valenti et al. 2015), Mg i λ1.53 µm is detected in the H-band range (Maguire et al. 2010b;Valenti et al. 2015;Yuan et al. 2016), while in the K-band range the Brackett γ is possibly blended with Na i (Dall'Ora et al. 2014). The low number and weakness of metal lines reduce the risk of systematics effects produced by differences in progenitor metallicity (e,g., Dessart et al. 2014;Anderson et al. 2016). Schmidt et al. (1992) had already pointed out the benefits of measuring distances to SNe II using near-IR photometry. However, at present, there have been very few systematic studies (e.g., Schmidt et al. 1992 for the EPM; Maguire et al. 2010a, de Jaeger et al. 2015 for the SCM). In particular, Maguire et al. (2010a) suggested that it may be possible to reduce the scatter in the Hubble Diagram (HD) to 0.1-0.15 mag (distance precision of 5-7 per cent) using near-IR instead of optical photometry. However, this result is based on the analysis of a set of 12 SNe II, 11 of them at z < 0.01, so being highly affected by peculiar velocities. To test this promising result, de Jaeger et al. (2015) applied the SCM to a set of 24 SNe II at 0.01 < z < 0.04, obtaining a HD rms of 0.28 mag (distance precision of 13 per cent) for the J-band and therefore questioning the improvements of the SCM distance precision using near-IR photometry.
The goal of this study is to investigate the PMM distance precision using near-IR photometry.
We organize our work as follows. In Section 2 we describe the photometric and spectroscopic data. In Section 3 we present the PMM developed in R14. In Section 4 we develop an algorithm to achieve nonparametric light curve fitting. In section 5 and 6 we compute Galactic total-toselective broadband extinction ratios and K-corrections for BVIJHK bands, respectively. In Section 7 we compute host galaxy total-to-selective broadband extinction ratios and host galaxy colour excesses through the analysis of colourcolour curves. In Section 8 we estimate expansion velocities and explosion epochs. In Section 9 we apply the PMM to our SN II sample, constructing HDs for BVIJH bands. Discussion about the PMM distance precision and systematics are in Section 10. In Section 11 we present our conclusions.

OBSERVATIONAL MATERIAL
We base our work on data obtained over the course of the Carnegie Type II Supernova Survey (CATS; PI: Hamuy, 2002Hamuy, -2003, a program whose main objective was to study nearby (z < 0.05) SNe II. Optical photometry and spectroscopy, along with some near-IR photometry, were obtained with the 1-m Swope, 2.5-m du Pont, and 6.5-m Magellan Baade and Clay telescopes at Las Campanas Observatory. A few additional optical images were obtained with the 0.9-m and 1.5-m telescopes at Cerro Tololo Inter-American Observatory. During the CATS survey, 34 SNe II were ob-served. Optical photometry and spectroscopy of these SNe II, along with the description of the data reduction, is presented in Galbany et al. (2016) and Gutiérrez et al. (2017), respectively. Next, we briefly summarize the general techniques used to obtain the near-IR photometric data, which will be released in a forthcoming publication.

Near-IR Photometric Data
The near-IR photometric observations were obtained with the JHK bands mounted in the Swope Telescope IR camera and the Wide Field IR Camera on the du Pont Telescope. Images where processed with a collection of IRAF 1 tasks. These include dark subtraction, flat field correction, sky subtraction, image registration and stacking. Instrumental magnitudes were obtained using the point spread function (PSF) technique, implemented in the SNOoPY 2 package. The near-IR magnitudes of the reference stars were calibrated using standard star fields obtained soon before or after the target field with an airmass similar to the target field.

Sample of Supernovae
Among the 34 SNe II observed over the course of the CATS survey, we select a subset of 10 SNe II which comply with the following requirements: (1) having at least two photometric measurements in the BVIJH bands at 35-75 d since explosion (see Section 9), and (2) having at least one measurement of the expansion velocity at an epoch covered by the photometry mentioned in point 1. To this sample, we add 14 SNe II from the literature. Table 1 lists our final sample of 24 SNe II, which includes the SN name (Column 1), the name of the host galaxy and its type (Column 2 and 3), the heliocentric SN redshift and its source (Column 4 and 5), host galaxy distance measured with Cepheids, the Tip of the Red Giant Branch (TRGB), or SN Ia (Column 6), Galactic colour excess (Column 7), and references for the data (Column 8). We also use optical and near-IR spectra of SNe II with the purpose of computing total-to-selective broadband extinction ratios (Section 5 and 7) and K-corrections (Section 6), and to estimate explosion epochs (Section 8.2).

PHOTOSPHERIC MAGNITUDE METHOD
The absolute magnitude of a SN II during the photospheric phase depends strongly on the temperature and the size of the photosphere (e.g., Kasen & Woosley 2009;R14;Pejcha & Prieto 2015). The latter can be estimated from the velocity of the material instantaneously at the photosphere (hereafter, photospheric velocity, v ph ) and the time since the SN explosion epoch t0, under the assumption of homologous expansion (e.g., Kirshner & Kwan 1974). R14 1 IRAF is distributed by the National Optical Astronomy Observatory, which is operated by the Association of Universities for Research in Astronomy (AURA) under cooperative agreement with the National Science Foundation. 2 SNOoPy is a package for SN photometry using PSF fitting and/or template subtraction developed by E. Cappellaro. A package description can be found at http://sngroup.oapd.inaf.it/snoopy.html found that the time since explosion works better than the V −I colour (used as a proxy for temperature) to standardize the brightness of SNe II (see Fig. 9 in R14), showing that for a given band x the absolute magnitude in any moment ti during the photospheric phase, Mx,∆t i ,v ph,i , can be parametrized as (1) Here, ∆ti ≡ (ti − t0)/(1 + z) is the elapsed time since the explosion in the SN rest frame at redshift z, and ax,∆t i is a function that can be calibrated empirically. Previously, Kasen & Woosley (2009) found similar results for the SN II brightness standardization, but using SN II models.
With the knowledge of t0 and a measurement of v ph in any stage of the photospheric phase, we can compute the absolute magnitude at ti (equation 1) and, therefore, compute the SN distance modulus given by m corr Here, mx,i is the apparent magnitude, A G x,i and A h x,i are the Galactic and host galaxy broadband extinction, respectively, and Kx,i is the K-correction. If more than one measurement of v ph is available, then we compute the distance modulus through a likelihood maximization (see Section 9.1).

LIGHT CURVE FITS
In equation (2) we need all quantities at the same epoch. Being more time consuming, spectroscopy is in general less abundant than photometry, so performing photometric interpolations is a reasonable choice. Previous efforts to fit SN II light curves use both parametric and nonparametric methods. Parametric methods assume parametric functions that capture the behaviour of the light curve from early to late stages, where parameters are obtained through least-square minimization (e.g., Olivares E. et al. 2010) or through Bayesian methodologies (e.g., Sanders et al. 2015). Nonparametric methods are based on nonparametric regressions (NPR) like local regressions (e.g., Olivares 2008) and Gaussian processes (e.g., de Jaeger et al. 2017a). Since the light curves of the SNe in our sample are in general well sampled, we prefer to use NPR methods for the light curve fitting, thus avoiding the use of heuristic models.
In this work we make use of loess, a NPR method that performs polynomial fits over local intervals along the domain (Cleveland et al. 1992). To perform a loess fit, we have to specify: (1) the class of the local polynomial, which can be linear or quadratic, (2) the smoothing parameter, which defines the neighbourhood size around each element of the independent variable, where data can be well approximated by the aforementioned local polynomial, and (3) the distribution of random errors, which can be normal or symmetric (for more details, see Cleveland et al. 1992). We assume the null hypothesis that residuals are normally distributed, which can be checked with a normality test. An optimal value for the smoothing parameter can be obtained from data using the "an" information criterion (AIC, Akaike 1974, see Appendix A). Therefore, to perform a loess fit, we only have to decide the local polynomial. We choose a  Leonard et al. (2002); S05: Sollerman et al. (2005); here: this work. ⋆ Host galaxy distance moduli measured with Cepheids (⊘), TRGB (⊗) with the Jang & Lee (2017b) calibration, or with SN Ia (△). ⋄ Galactic colour excesses from Schlafly & Finkbeiner (2011), with an error of 16 per cent (Schlegel et al. 1998). ⊳ Saha et al. (2006)  quadratic polynomial in order to give more freedom to the loess fitting procedure. When the loess routine cannot perform a fit (e.g., for light curves with less than six points), we perform a low-order (linear or quadratic) polynomial fit.
To test whether photometry errors can account for the observed dispersion around the light curve fit, f fit x,t (in flux units), we compute its log-likelihood given by where fx,i and σ f x,i are the apparent magnitude and its error in flux units, respectively, and σx,int is the intrinsic error. If an intrinsic error is necessary to maximize the loglikelihood, then we add it in quadrature to the photometry errors and perform again the light curve fitting. We repeat this process until an intrinsic error is not necessary.
To test the normality of the residuals, we use the Rescaled Moment (RM; Imon 2003) test (see Appendix A). Among all light curves fits, 80 per cent have residuals with RM p-values ≥ 0.05, for which the hypothesis that residual are normally distributed cannot be rejected within a confidence level (CL) of 95 per cent. For the remaining 20 per cent, light curve fits are still unbiased and consistent, but the confidence interval (CI) of the parameters may be un-trustworthy (Doane & Seward 2016). Anyway, in this work, to prevent any shortcoming related to the non-normality of the residuals, we perform simulations to compute CIs.
To compute the CI around a light curve fit, we perform 10 4 simulations varying randomly the photometry according to its error. For each realization we perform a loess (or a low-order polynomial) fit, thus obtaining 10 4 simulated light curves per band. These simulations will allow us to compute its probability density function (pdf) at different epochs. Fig. 1 shows results of the aforementioned fitting procedure applied to the SN II light curves used in this work, where solid lines are the loess (or low-order polynomial) fits, while shaded regions indicate values between the 10th and the 90th percentile, i.e., the 80 per cent CI.

GALACTIC BROADBAND EXTINCTION
In equation (3), the Galactic broadband extinction in a photometric band x is given by  (Olivares 2008). Here, λo is the wavelength in the observer's frame, λr = λo/(1 + z) is the wavelength in the SN rest frame at redshift z, S x,λo is the x-band transmission function, F λr,i is the spectral energy distribution (SED) of the SN at epoch ti. A G λo and A h λr are the Galactic and host galaxy monochromatic extinctions, respectively, given by where R λ ≡ A λ /E(B −V ) is the extinction curve for our Galaxy (R G λ ) and hosts (R h λ ), and EG(B −V ) and E h (B −V ) are the Galactic and host galaxy colour excess, respectively.
Since the SED of SNe II evolve with time, we expect that the broadband extinction A G x,i also evolve with time 3 . As the SED of SNe II has a blackbody nature, hereafter we 3 We remark the difference between a monochromatic extinction A λ , which is constant for a fixed wavelength λ, and a broadband use the intrinsic B −V colour (a proxy for temperature) to represent its evolution.
In a previous work, Olivares (2008) In this work, in order to convert Galactic colour excesses directly into Galactic broadband extinctions suitable for SNe II, we compute the Galactic totalto-selective broadband extinction ratios R G x,i , such that With the purpose of obtaining representative R G x,i values for a local SN II sample through equation (8) and (5), we use: (1) a library of dereddened and deredshifted SN II spectra (see Appendix B), (2) colour excesses and redshifts from the following representative ranges: EG(B −V ) = 0.0-0.36 mag, E h (B −V ) = 0.0-0.83 mag, which were taken from the extinction A x,i , which depends on the SED and the x-band transmission function.  Figure 2. Galactic total-to-selective broadband extinction ratios for BVI (left) and JHK (right) for SNe II as a function of the intrinsic B −V colour, along with residuals. Solid black lines correspond to the polynomial fits, red short-dashed lines correspond to loess regressions to the residuals, while blue long-dashed lines correspond to residuals between the m15mlt3 model (Dessart et al. 2013) and polynomial fits from observations. Gray regions indicate values within one rms, while black dash-dotted lines are the inner fences.
SN sample reported in R14, and z = 0.0-0.032, and (3) an extinction curve to redden the spectra for both our Galaxy and hosts. For the latter, since a representative extinction curve along the line of sight of SNe II is still controversial, we adopt the Fitzpatrick (1999) extinction curve with RV = 3.1. For each spectrum, we perform 10 4 simulations picking randomly values of EG(B −V ), E h (B −V ), and z from the aforementioned ranges, adopting the median as the R G x,i representative value and the 80 per cent CI as its error. The left side of Fig. 2 shows the R G x,i values as a function of B −V for BVI bands. There is a clear dependence of The y-axis scale at the left of Fig. 2 is the same in the three panels, so we can see that the redder the band the less the dependence on B −V , with R G I,i being nearly constant. This behaviour is due to the blackbody nature of the SN II SED, where for the longer wavelengths the less the dependence of the SED slope on temperature.
To express the dependence of R G x,i on B −V we perform polynomial fits. The latter, unlike NPR methods like loess, allow to perform corrections in an easy and less timeconsuming way (see Appendix C).
To determine the optimal degree for the polynomial fit, we consider two criteria: the AIC and the Bayesian information criterion (BIC, Schwarz 1978). For more details, see Appendix A. Based on evidence ratios (Table F1), for the B-, V -, and I-band, the AIC favours degrees ≥ 3, ≥ 1, and ≥ 1, respectively, while the BIC favours degrees between 2 and 6, 1 and 3, and 1 and 4, respectively. Results for both criteria are consistent. By the principle of parsimony (a.k.a. Occam's razor), we adopt the lowest degrees, i.e., 2, 1, and 1 for the B-, V -, and I-band, respectively. For JHK bands (right of Fig. 2) we adopt constant values. Although the small number of near-IR spectra means that the results are not fully statistically robust, we are confident about the negligible dependence of R G x,i on B −V for JHK bands based on the small rms values we obtained ( 0.001).
Once the optimal polynomial degrees for R G x versus B −V are determined, we perform 10 4 bootstrap resampling of the data in order to compute the polynomial fit parameters and their errors, adopting the median as the representative value. Results are summarized in Table F2.
The bottom of each panel in Fig. 2 shows the residuals of the polynomial fit. To identify possible outliers we use the Tukey (1977) rule, where values below Q1−1.5(Q3−Q1) or above Q3 + 1.5(Q3 − Q1) (known as inner fences, where Q1 and Q3 are the first and third quartile, respectively) are considered outliers. The few points detected as outliers are consistent with being within inner fences considering their errors, so we do not discard them from the analysis. To analyse possible trends not captured by the polynomial fit, we perform a loess regression (red short-dashed line) to the residuals. Variations in the loess fits are mostly within one rms, meaning that the evolution of R G x on B −V can be well represented by a polynomial fit of degree determined with the AIC/BIC. For all bands we obtain RM p-values > 0.05, which means that we cannot reject the null hypothesis that residual are normally distributed (95 per cent CL). Based on this, we can treat the R G x rms error as a normal one.
For comparison, we compute R G x for BVIJHK bands using synthetic spectra of the m15mlt3 model of Dessart et al. (2013). Residuals between the m15mlt3 model and polyno-mial fits from observations (blue long-dashed lines in Fig. 2) are mostly contained within one rms.

K-CORRECTION
The K-correction in a photometric band x is given by (Olivares 2008), being K s x,i the selective term. We proceed in the same way than in Section 5, but now the evolving SED is modified by SN host galaxy colour excess and redshift.
As in Section 5 we aim for an analytical expression for K s x , for which we perform polynomial surface fits as a function of B −V and z. Since Kx = 0 for z = 0, any zindependent term on the K s x polynomial fit is zero. Dividing by z, the polynomial surface to adjust will be of the form being OB−V and Oz,j 1 the orders in B −V and z, respectively, and aj 1 ,j 2 the fit parameters.
To determine the orders OB−V and Oz,i, we generate 10 5 spectral samples, where for each sample we assign to each spectrum a random redshift up to 0.032. For each realization, we obtain optimal order values using the AIC/BIC and the principle of parsimony. In all cases we obtain that K s x /z depends only on B −V , i.e., it is z-independent for z ≤ 0.032. Fig. 3 shows the K s x /z values as a function of B −V for BVI (left) and JHK (right). We perform the same analysis than in Section 5. For BVI bands we adopt straight lines, while for JHK bands we fit constant values (see Table F1). Results are summarized in Table F2. Variations of the loess fits to the residuals are within one rms, meaning that the dependence of K s x /z on B −V can be well represented by the polynomial fit of degree determined with the AIC/BIC. For all bands we obtain RM p-values > 0.05, which means that we cannot reject the null hypothesis that residual are normally distributed (95 per cent CL). Based on the latter, we can treat the K s x /z rms error as a normal one.

HOST GALAXY BROADBAND EXTINCTION
In equation (3), the host galaxy broadband extinction in a photometric band x is given by (Olivares 2008). We proceed in the same way than in Section 5 and 6, but now the evolving SED is modified only by the SN host galaxy colour excess. Similar to Section 5, we define the host galaxy total-toselective broadband extinction ratios R h x,i , such that The optimum R h x versus B −V polynomials and their parameters are summarized in Table F2.
R14 showed that for SNe II the B −V versus V −I colour-colour curve (C3) can be used to estimate E h (B −V ) through the method proposed by Natali et al. (1994), which was originally developed to estimate interstellar colour excess for open clusters.
The C3 method states that, under the assumptions that (1) the C3 can be well-represented by a straight line, and (2) all SNe II have the same C3 (which means the same slope and intercept), the host galaxy colour excess can be estimated with the formula (e.g., Munari & Carraro 1996;R14). Here, S ≡ {cx, cy} indicates the colours used as x-and y-axis in the colourcolour diagram, corrected for Galactic colour excess and K-correction. M{mS} is the median of a set of SN II C3 slopes {mS}, nS,i = cy,i − M{mS} × cx,i and nS,zp are the y-intercept of the C3 linear fit using a fixed slope M{mS} and that of the SN II less affected by colour excess, respec- (14) must be evaluated separately at each point of the C3. In principle, one colourcolour observation is enough to estimate the colour excess, however more observations allow to check internal consistency and reduce observational errors.
The C3 method relies strongly on the aforementioned two assumptions. In a previous work, R14 assumed the linearity of C3s based on the blackbody nature of the SED of SNe II during the photospheric phase, while assumption 2 was adopted based on the dependence of the emergent flux mainly on temperature displayed by SN II atmosphere models (e.g., Eastman et al. 1996;Jones et al. 2009). In this work we show that C3s can indeed be expressed as straight lines for several colour combinations (see Appendix C). Therefore, the major source of systematics comes from assumption 2. There are indeed some effects, like the line blanketing, that modify the SED continuum shape. In addition, differences in photometric systems (Scorrection; Stritzinger et al. 2002) are expected to produce further changes on C3 parameters. Therefore it is propitious to search for a colour combination where the effect of a colour excess on a C3 is greater than the effect of systematics.
An analysis of the effect of systematics on the C3 yintercept is beyond the scope of this work because it requires an ample set of unreddened SNe II. However, the effect of systematics on C3 slopes and its consequent effects on the E h (B −V ) estimation through the C3 method can be quantified in a simple way.
The presence of dust along the line of sight produces a vertical displacement of the C3 (for a graphical representation, see R14) where, following equation (14), the magnitude of the displacement and its rms error are respectively. In equation (16), we do not include the error . Selective term of the K-correction over redshift for BVI (left) and JHK (right) for SNe II as a function of the intrinsic B −V colour, along with residuals. Solid black lines correspond to the polynomial fits. Red short-dashed lines correspond to loess regressions to the residuals, while blue long-dashed lines correspond to residuals between the m15mlt3 model (Dessart et al. 2013) and polynomial fits from observations. Gray regions indicate values within one rms, while black dash-dotted lines are the inner fences.
induced by errors in γS,i, which is lower that 17 per cent of the uncertainty induced by the error in M{mS}. In order to find the colour combination that maximizes the dust effect (equation 15) and minimizes its error (equation 16), we define the quantity (a signal-to-noise ratio) Therefore, the most appropriate colour set S to compute E h (B −V ) with the C3 method is one that maximizes ξS,i. Fig. 4 shows the ξS,t values for all possible independent colours combinations with the BVIJH bands, using the M{mS} and rms{mS} values computed with our SN II sample (see Appendix C), and using B −V = 0.0 and 1.4, which are typical colours at the start and end of the photospheric phase, respectively. We do not include the K-band in this analysis because of the scarcity of SNe II with photometry in that band. The best colours combinations, independent of the intrinsic B −V , are those involving B −V , with V −I versus B −V the best. For this combination we obtain M{mS} = 0.45 ± 0.07. We remark that colours combinations that do not include the B-band have ξS,t 1.0, which indicates that the noise induced by intrinsic differences of C3 slopes is greater than the effect of host galaxy dust, and therefore those combinations are not suitable for E h (B −V ) measurement through the C3 method. We point out that colours combinations under the diagonal correspond to those above the diagonal but with the axes exchange. In principle they give the same information. However, by construction, they maximize displacement in x-axis instead of y-axis.
To compute the pdf of nS,zp for S = {B −V , V −I}, we use the data of SN 2003bn and SN 2013ej, which are affected  by a negligible host galaxy colour excess (R14), maximizing the likelihood of a straight line with slope 0.45 ± 0.07. With this process, we obtain a pdf with median of 0.107 mag and rms = 0.053 mag. Since the RM p-value for the latter distribution is > 0.05, we treat it as a normal distribution.
#SNe=13 offset=-0.01 mag rms=0.05 mag  Table F3), which are based on the fit between observed spectra and SN II models. We measure a median offset of −0.01 mag, meaning that our estimations of E h (B −V ) are slightly lower than those estimated by Olivares E. et al. (2010). Both methods are consistent within ±0.05 mag.

EXPLOSION EPOCH AND PHOTOSPHERIC VELOCITY
The explosion epoch and the photospheric velocity are, under the assumption of homologous expansion, the unique parameters determining the actual size of the photosphere (Kirshner & Kwan 1974).

Photospheric Velocities
The most widely used method to estimate SN photospheric velocities consists of measuring the blueshift of P Cygni absorption minima in SN spectra (Kirshner & Kwan 1974;Eastman & Kirshner 1989). Weak lines, like those from Fe ii species, are typically used under the assumption that they are formed near the photosphere (e.g., Leonard et al. 2002). A more confident method to estimate photospheric velocities is through the cross-correlation technique Takáts & Vinkó 2012), where observed spectra are compared to those from SN models which have known photospheric velocities. The application of the latter method is beyond the scope of this work, therefore we will use velocities derived from the Fe ii λ5169 line absorption minima as a proxy for the photospheric velocity.
To estimate Fe ii λ5169 absorption minima with appropriate errors, we have to consider the uncertainties induced by the noise and spectral resolution (∆λ) of each spectrum, and also by the endpoints we choose for the line profile.
We estimate the noise on the Fe ii λ5169 line profile of each spectrum performing a loess fit and then removing it to the observed line profile. Then we generate 10 4 simulated line profiles, varying randomly the noise over the loess fit, wavelengths within ∆λ, and endpoints. For each realization we apply a loess fit, registering the minimum value. The output of this process is a distribution of absorption minima, which we convert to velocities using the relativistic Doppler equation. With this process we obtain typical v ph rms errors between 30 and 230 km s −1 , with a median of 76 km s −1 .
Photospheric velocities are estimated from spectroscopic data, corrected for the SN heliocentric redshift.
In some cases, SN II spectra shows narrow emission lines as result of a superposed H ii region at the SN position. These narrow lines allow a good estimation of the SN heliocentric redshift, under the assumption that the SN is spatially close to the H ii region (e.g., Anderson et al. 2014a). When those lines are not present in the SN spectra, the heliocentric redshift of the host galaxy is used as a proxy for the SN heliocentric redshift. However, since most of the SNe II in our set explode in spiral galaxies, the SN heliocentric redshift has a component due to the galaxy rotation. Anderson et al. (2014a) computed heliocentric redshifts of 72 SNe II using H ii region narrow emission lines, and comparing with heliocentric redshifts of the host galaxy nucleus, they obtained a zero-centred distribution with a rms of 162 km s −1 , which is attributed to the galaxy rotation effect.
In our sample, 11 SNe II (SN 2002gd, SN 2002gw, SN 2002hj, SN 2003B, SN 2003E, SN 2003bl, SN 2003bn, SN 2003ci, SN 2005ay, SN 2009N, and SN 2014G) show H ii region narrow emission lines in the spectra, which we use to estimate the heliocentric redshift. Another six SNe (SN 2003T, SN 2004et, SN 2005cs, SN 2008in, SN 2012aw, and SN 2013ej) exploded within nearly face-on galaxies, in which case we adopt the redshift of the host galaxy nucleus. For SN 1999em we adopt the value from Leonard et al. (2002), and for SN 2003hn we use the average of the Na i D velocities measured by Sollerman et al. (2005). The remaining five SNe (SN 2003cn, SN 2009ib, SN 2009md, SN 2012A, and SN 2012ec) did not occur within nearly face-on galaxies, and do not show H ii region narrow emission lines in the spectra. For those cases we adopt the redshift of the host galaxy nucleus, with an error of 162 km s −1 (that we assume normal) to take into account the host galaxy rotational velocity. Adopted SN heliocentric redshifts are listed in Table 1.

Explosion Epoch
The SN explosion epoch can be estimated by means of photometric information; it can be constrained between the last non-detection t ln and the first detection t fd (e.g., Nugent Valenti et al. 2016), or estimated through a polynomial fit to the rise-time photometry when it is available (e.g., González-Gaitán et al. 2015;. The spectroscopy of a SN can also provide information about its explosion epoch by means of the comparison with other spectra of SNe with explosion epoch estimated through photometric information (e.g., Anderson et al. 2014b;Gutiérrez et al. 2017).
Column 4 and 5 of Table F3 lists the t ln and t fd values of the SNe in our set, respectively. The explosion epochs for our set are typically constrained within 14 d using photometric information, which is twice the range suggested by R14 (namely, 7 d) to reduce errors induced by t0 errors over PMM distances. We need, therefore, to include spectroscopic information in order to better constrain the explosion epochs.
As was done by Anderson et al. (2014b) and Gutiérrez et al. (2017), to estimate t0 we use optical spectroscopy along with the Supernova Identificator code (SNID; Blondin & Tonry 2007), which finds by crosscorrelation the spectra from its SN library that are more similar to the input spectrum. For a good estimation of t0 with SNID, we need a library with spectra of an ample amount of SNe II that sample the high spectral diversity displayed by SNe II (e.g., Gutiérrez et al. 2017) and with t0 constrained by photometric information. In this work, we compile optical spectroscopy of 59 SNe II with t0 constrained within 10 d (for more details, see Appendix D).
To estimate the explosion epoch of a given SN (SNinput) with N spectra ({spec}) through SNID and using our SN II templates library, we perform the following procedure: 1. We run SNID using as input the N spectra of SNinput earlier than 40 d since the first detection. The SNID output for each spectrum is a list with the best-matching templates, their phase since explosion, and their rlap parameter (which indicates the strength of the correlation).
2. We convert phases since explosion to explosion epochs (since we know the phase of each SNinput spectrum). The associated errors are derived from the rlap values through a procedure described in Appendix D.
3. From each of the N lists, we select the first ten bestmatching templates with rlap > 5.0, compiling them in a unique list. From this list, we extract a sublist for each of the M best-matching SNe (SN bm ). With each of the M sublist, we compute the SNin explosion epoch as the average, taking the standard deviation as the associated error, and including the respective explosion epoch error of the SN bm through a Monte Carlo error propagation. If a spectrum gives a median t0 greater than 40 d, then we remove it from the analysis.
4. After that, we compute the likelihood L(t0|{spec}) with the M results, including an error of 4.1 d which is the rms obtained from the comparison between explosion epochs constrained with photometric information and those derived with SNID (see Appendix D).

APPLYING THE PMM
Once all observables required for the PMM are available, the next step is to prepare the data before applying the method. As was mentioned in Section 4, we interpolate photometry to the epochs of the photospheric velocities. Since we want to study the PMM distance precision at different photometric bands, i.e., changing only the photometry, we use epochs where spectroscopy is covered simultaneously by optical and near-IR photometry. In the case of SN 2002hj, it does not have spectroscopy covered by J-band photometry, so we interpolate photospheric velocities (using loess) and BVIH photometry to the epochs of J-band photometry.

Calibration
For the PMM calibration, we express ax,∆t as where ZPx is the zero-point of the PMM in the x-band, and a * x,∆t is a function that represents the dependence of ax,∆t on ∆t (without the constant term).
To estimate the evolution of a * x,∆t with ∆t, we use the a * x,∆t i values of the SNe in our set with two or more v ph measurements during the photospheric phase. For each SN, the a * x,∆t i values are given by where δSN is an additive term to normalize the a * x,∆t i values of each SN to the same scale. Based on the definition given in R14, the dependence of a * x,∆t on ∆t has the form a * x,∆t = fx,∆t − 5 log ∆t 100 d .
We express the dependence of fx,∆t on ∆t through polynomials. We use the AIC/BIC to determine the optimum polynomial order for fx,∆t and the values of δSN, while to estimate the time range of applicability of the PMM, we group the fx,∆t i values in bins of width 10 d and then we compute the rms of the points in each bin. We found that rms values are lower in a range 35-75 d since the explosion. Among all optimum orders for BVIJH bands (see Table F1), we select the order that the different bands have in common, i.e., order one. With this, we prevent that differences in the rms(fx,∆t) value for different bands are due to differences in the order of the polynomial fit. To estimate error in the parameters, we perform 10 4 bootstrap resampling. Table 2 lists fx,∆t fits parameters for BVIJH bands.
The left half of Fig. 6 shows the values of a * x,∆t i as a function of ∆ti for BVIJH bands. The variation of the loess fits (red dashed lines) are within one rms (black dotted lines), which means that polynomial fits we adopted capture almost all the dependence on ∆t.
The PMM zero-points can be obtained using a sample of SNe II at known distances where, for each SN, we have Here, µ SN host is the SN host galaxy distance modulus and µ * x is the SN pseudo-distance modulus. The latter, for each measurement of v ph at time ti, is defined similar to equation 2 but with a * x,∆t i instead of ax,∆t i , i.e., The pdfs of µ * x,i are obtained through equation 23 using the pdfs of the observables for each photospheric velocity epoch. Finally, we combine the pdfs of µ * x,i in a unique µ * x pdf maximizing the likelihood (equation A1) for a constantonly model. To compute accurate ZPx values, we need SNe II in galaxies with distances measured with the best possible precision. Among the SNe that we compiled from the literature, there are only three (SN 1999em, SN 2003hn, and SN 2012aw) in galaxies with distances measured through Cepheid, and five (SN 2003hn, SN 2004et, SN 2005cs, SN 2012aw, and SN 2013ej) in galaxies with distances measured with TRGB. Cepheid distances for the hosts of SN 1999em and SN 2012aw were reported by Saha et al. (2006), while Riess et al. (2016) reported the Cepheid distance of the host of SN 2003hn. Comparing Cepheid distances of six galaxies in common between the two publications (NGC 1365, NGC 3370, NGC 3982, NGC 4536, NGC 4639, and NGC 5457) we found that Cepheid distances reported by Saha et al. (2006) are, on average, 0.19 mag greater than those reported by Riess et al. (2016), showing a rms of 0.13 mag. The latter could indicate a systematic difference between the two calibrations, which can introduce an undesirable noise on the ZPx estimation if we rescale Saha et al. (2006) distances to the Riess et al. (2016) calibration. For this reason, we decide to use only SNe in galaxies with TRGB distances, that can be homogenized to the Jang & Lee (2017a) calibration, which is based on the distance to the Large Magellanic Cloud and NGC 4258. Recalibrated TRGB distances are listed in Column 6 of Table 1. From these five SNe II, we discard SN 2004et since the TRGB distance of its host is at least 0.59 mag higher compared to the distances we compute for SN 2004et and two other SNe II that exploded in the same galaxy (see Appendix E).
The right half of Fig. 6 shows the ZP SN x values for BVIJH. As in the case of µ * x , the pdf of ZPx is obtained combining the pdfs of ZP SN x . Median values, 99 per cent CI, and rms values for ZPx are summarized in Table 2.
Once the PMM zero-points are computed, we can estimate the distance modulus for each band as µx = µ * x − ZPx. Median values, 80 per cent CI, and rms values for µx are summarized in Table F4, where we include the TRGB zeropoint systematic error of 0.058 mag (Jang & Lee 2017a).

Hubble Diagrams
To investigate the PMM distance precision, we construct HDs. We convert heliocentric host galaxy redshifts to cosmo- logical ones using as reference the cosmic microwave background (CMB) dipole (Fixsen et al. 1996). Redshift errors are dominated by peculiar velocities, with a rms of 382 km s −1 for local SN Ia host galaxies (z < 0.08, Wang et al. 2006), followed by the error in the determination of the Local Group velocity (rms of 187 km s −1 , Tonry et al. 2000). CMB redshifts and their rms errors are listed in Table F4.
Taking into account the pdfs of the pseudo-distance moduli (µ * x ) and the pdfs of the CMB redshifts (czCMB), we compute the Hubble diagram intercept (aHD,x) maximizing the likelihood (equation A1), where the model for the pseudo-distance modulus is given by the Hubble law µ * x = aHD,x + 5 log(czCMB).
The left half of Fig. 7 shows HDs for BVIJH bands, using PMM distances for all SNe in our set. The rms, greater than 0.5 mag for all bands, is mostly produce by peculiar velocities of host galaxies at low redshift. In fact, the median redshift of the host galaxies in the HD is 1528 km s −1 , where a redshift error of 382 km s −1 translates into a magnitude error of 0.54 mag. Indeed, if we use redshifts corrected for the infall of the Local Group toward the Virgo cluster and the Great Attractor (czcorr) instead of CMB redshifts, we obtain a HD rms of 0.34-0.38 mag for VIJH bands (see the right half of Fig. 7). We note that even after infall corrections the scatter in the HDs is still mostly due to SNe in galaxies with cz < 2000 km s −1 . Therefore, to estimate the PMM distance precision and the Hubble constant (H0), given by log H0 = (25 − aHD,x + ZPx)/5, we use only SNe II with cz > 2000 km s −1 and, as visible in the left half of Fig. 8, the HD rms decreases significantly. The corresponding values of H0 and rms are listed in Table 3. The values of H0 range between 67.1 and 74.9 km s −1 Mpc −1 . Taking into account that the ZPx values were calibrated using TRGB distances in the scale of Jang & Lee (2017a), our H0 values, as expected, are consistent within the errors with those reported in Jang & Lee (2017b), i.e., 71.17 ± 1.66 ± 1.87 km s −1 Mpc −1 , which also use the Jang & Lee (2017a) calibration.
As visible in Column 4 of Table 3, all the H0 values are compatible within their errors. However, we note that H0 decreases moving from shorter to longer wavelengths, which  could suggest a systematic introduced by: (1) our assumption of the RV value for the SN host galaxies, or (2) an underestimation and/or an overestimation of the host galaxy colour excesses for the four SNe in the PMM calibration set and the the nine SNe at czCMB > 2000 km s −1 , respectively. To test the first possibility, we change the RV adopted for SN II host galaxies to the lowest RV value for which the Fitzpatrick (1999) extinction curve is defined (RV = 2.3). As is visible in Column 6 of Table 3, there are no significant changes in the H0 values. For the second possibility, we found that an underestimation of E h (B −V ) for the SNe in the calibration set, or an overestimation of E h (B −V ) for the SNe at czCMB > 2000 km s −1 , of 0.05-0.07 mag can erase the tension between the H0 values for all bands. Given the typical statistical E h (B −V ) error of 0.097 mag (see Section 7), the probability of obtaining an E h (B −V ) underestimation of 0.05-0.07 mag with four objects is of 8 per cent, while to obtain an overestimation in a same amount for nine objects is of 5 per cent. It is worth mentioning that also the scatter in ZPx decreases going from shorter to longer wavelengths, suggesting again a trend introduced by the combination of a large uncertainty in the colour excess estimation and the low number statistics.
Regarding the HD scatter, we note that the rms of 0.28 mag obtained with the B-band decreases to 0.15-0.18 mag for the VIJH bands. Despite the good results, the low number of SNe II within galaxies at czCMB > 2000 km s −1 makes the result not statistically robust. Therefore, to check the PMM distance precision, it is necessary to include more SNe II into the analysis. Thus, we include the four SNe II used for the PMM calibration, plus other two in galaxies having Cepheid and SN Ia distances. The latter distances are converted to redshifts through the Hubble law (equation 24), where for each band we use the H0 value listed in Column 4 of Table 3 and the ZP value given in Table 2.
The right half of Fig. 7 shows the HDs computed with the selected SNe II based on the aforementioned criterion, which correspond to our final sample. For VI bands we obtain a HD rms of 0.15-0.16 mag. The lowest HD rms is obtained with the J-band, whose rms of 0.13 mag translates into a distance precision of 6 per cent. This value, compared to the rms of 0.15-0.26 mag obtained for optical bands, suggests that using the J-band photometry instead of optical one to measure PMM distances can improve the precision by at least 0.07 mag.
For the H-band we expected a similar HD rms than for J-band since, among BVIJH bands, the H-band is the least affected by dust extinction. We, however, obtained a HD rms of 0.15 mag. The latter can be partially due to the higher photometry error of the H-band (of 0.07 mag) with respect to the error of the J-band (of 0.05 mag).

Statistical significance of the result
Given the small size of our SN sample, the HD rms of 0.13 mag we measured for the J-band is not statistically robust enough to be considered as a measure of the PMM distance precision in that band. In particular, we want to know the probability of measuring a rms ≤ 0.13 mag with N = 15 values drawn from a parent distribution with standard deviation (σ) ≥ 0.13 mag. Assuming a normal parent distribution, the quantity (N −1)(rms/σ) 2 has a chi-square distribution with N − 1 degrees of freedom. Using this, we found that there is a 1 per cent chance that the parent distribution has σ = 0.23 mag. Therefore, with the evidence we have, we can set an upper limit on the PMM distance precision with the J-band of 10 per cent with a CL of 99 per cent.

Comparison with other methods
For a consistent comparison of our results with those from other SN II distance measurement methods, we select results from works that uses a sample of SNe II at z ≈ 0.01-0.03. Table 4 lists the best distance precisions reached by other SN II distance measurement methods along with results obtained in this work. We note that the precision we report in this work is significantly lower than the best dispersion obtained by other works with SCM and PCM.
We also compare PMM and SCM applied to the same sample for J-and H-band. For this comparison we discard SN 2002hj because there is not photometry in J-band at 50 d since explosion. As visible in Table 4, the dispersion is lower by ∼30-40 per cent in both band.

Error budget
Taking into account that the lowest HD rms is obtained with the J-band, in Table 5 we show the statistical error budget for the distances measured for that band. We see that 88.6 per cent of the statistical error is induced by the errors of the first four terms: the host galaxy colour excess, the ex-  plosion epoch, the PMM zero-point, and the J-band photometry. Therefore, it is possible improve the performance of the PMM in the future developing a better method to estimate E h (B −V ), selecting SNe II with explosion epoch constrained within a small range of time, including more SNe II in the PMM calibration set, and increasing the quality of the J-band photometry.

Diminishing systematics
Our results show that we are reaching a precision in distance modulus of ±0.1 mag with the PMM at near-IR wavelengths. Therefore it is important to control systematics, and push them below 0.1 mag. For the latter in the following, we analyse sources of systematics affecting our results: Fig. 6) is stronger at early times, so t0 errors have a strong effect at those epochs. In order to obtain errors lower than 0.1 mag for observations at ∆t 35 d, we need SNe II with explosion epochs constrained within 5 d.

Explosion epoch: The dependence of the PMM calibration on explosion epoch (left half of
2. SN heliocentric redshift: When the host galaxy heliocentric redshift is used as a proxy of the SN heliocentric redshift, a systematic error of is introduced into the absolute magnitude (equation 1). This effect increases when the photospheric velocity decreases, translating into errors 0.1 mag for photospheric velocities 3500 km s −1 . Therefore, if optical spectra of a SN II do not show H ii narrow emission lines due to a nearby H ii region or if the SN is not within a nearly face-on galaxy, then epochs for which photospheric velocities are greater than 3500 km s −1 are preferable.
3. Host galaxy redshift: A galaxy is believed to be within the Hubble flow when its redshift is greater than 0.01. At that redshift, peculiar velocities are thought to be negligible compared with the velocities due to the Universe expansion. However, the typical error of 382 km s −1 translates into a distance modulus error of 0.28 mag for z = 0.01. Including the error in the determination of the Local Group velocity of 187 km s −1 , the redshift error increases to 425 km s −1 . Therefore, in order to reduce the error induced by redshift errors at a level lower than 0.1 mag, it is necessary to observe SNe II within galaxies at z > 0.03.

CONCLUSIONS
Our main results are the following: 1. Using nine SNe II at cz > 2000 km s −1 , we obtained H0 ranging between 67.1 and 74.9 km s −1 Mpc −1 , and a HD rms of 0.15-0.28 mag.
2. Adding six SNe II with host galaxy distances measured with TRGB, Cepheids, or SN Ia (total 15), which distances were converted to redshifts through the Hubble law, we obtain a HD rms of 0.15-0.26 mag in the optical bands, which reduces to 0.13 mag in the J-band.
In order to test further the promising results we are obtaining in this work, it is necessary to carry out an optical and near-IR photometric follow-up of SNe II at z > 0.03 and with explosion epochs constrained within 5 d. For these SNe, it is necessary to take at least one optical spectrum at any epoch between 35-75 d since explosion.
Its evident from Fig. 1 that the quality of the near-IR photometry used in this work is in general lower than the optical one. Therefore, we expect that increasing the quality of the near-IR photometry will further improve our results.

APPENDIX A: MODEL SELECTION
For the model selection we consider two criteria: the "an information criterion" (AIC, Akaike 1974), which is based on information theory, and the Bayesian information criterion (BIC, Schwarz 1978), which is based on Bayesian inference. From a set of R models, the AIC selects the model that have the least information loss with respect to the unknown true model, while the BIC selects the model with the highest likelihood L, given by Here, Xi is the i-th observed data, p(Xi|model) is the probability density function (pdf) of Xi given the model, and N is the number of observed data points. Quantitatively, the AIC and BIC search for a balance between overfitting and underfitting penalizing the likelihood. For the AIC and the BIC, the best model is one which minimizes the quantity (corrected for small sample sizes, Sugiura 1978), and BIC ≡ −2 ln Lmax + k ln N (A3) (Schwarz 1978), respectively, where Lmax is the maximum likelihood achievable by the model, and k is the number of parameters of the model. It is known that a model selection based only on the minimum AIC value reached for a certain model does not provide enough evidence to prefer that model over the other ones (e.g., Akaike 1978;Burnham & Anderson 2002). Instead, it is necessary to include into the analysis the strength of evidence in favour of each model. To quantify the latter, it has been proposed to use the likelihood of the model given the data (e.g., Akaike 1978) which, normalized by the sum of likelihoods of all R models, defines the Akaike weights (e.g., Burnham & Anderson 2002). Here, AICmin is the minimum AIC value reached among the R models used in the analysis. The same idea is applicable for the BIC (Burnham & Anderson 2002), which defines the Bayesian weights For the AIC and BIC, the evidence ratios defined as wi/wj and pi/pj, respectively, allow comparison of the evidence in favour of the i-th model as the best model versus the j-th model. As reference, if evidence ratios are greater than 13.0, then there is a strong evidence in favour of the i-th model over the j-th model (e.g., Liddle 2007). If several models have substantial support as the best (e.g., evidence ratios < 13.0), then, by the principle of parsimony, we select the one with less parameters.
In the case of least-square regressions, with random errors normally distributed and with constant variance, −2 ln Lmax = N ln 2πσ 2 + N (A6) (e.g., Burnham & Anderson 2002, p. 17), whereσ 2 is the average of squared residuals around the model. The AIC and BIC in this case can be expressed as Sinceσ 2 is computed from data, it must be considered as a model parameter.
In the case of nonparametric regressions (NPR), like loess, k is not defined. Instead, it has been proposed to use the trace of the smoother matrix, tr(H), which is a quantity analogous to the number of parameters in a parametric regression (Cleveland et al. 1992;Hurvich et al. 1998). Replacing k = tr(H) + 1 in equation (A7) allows us to obtain the AIC version for NPR presented by Hurvich et al. (1998).
To check the normality of random errors, it is necessary to carry out a normality test. As we do not measure random errors directly, we use residuals instead. However, widely used normality tests like the Shapiro & Wilk (1965) and the Jarque & Bera (1987) test, when applied over residuals, have little power to reject the null hypothesis even when the random errors are not normal (Das & Imon 2016). Imon (2003) proposed a statistic more suitable to verify normality for regressions, which is based on the Jarque & Bera (1987) test and on the use of unbiased moments. The statistic of the test, called Rescaled Moment (RM), is defined as RM ≡ N c 3 μ 2 3 /μ 3 2 /6 + c · (μ4/μ 2 2 − 3) 2 /24 (A9) (Rana et al. 2009), whereμi is the i-th sample moment, and c ≡ N/(N − k). Under the null hypothesis of a normal distribution, the RM statistic is asymptotically distributed as a chi-square distribution with two degrees of freedom.

APPENDIX B: SN II SPECTRA LIBRARY
In order to compute total-to-selective broadband extinction ratios (Section 5 and 7) and K-corrections (Section 6) for SNe II, it is necessary to know their SED. The latter can be estimated through theoretical models (e.g., Sanders et al. 2015;de Jaeger et al. 2015de Jaeger et al. , 2017b or, as in Olivares (2008) and in this work, through spectroscopic observations. In practice, spectra are not always taken with the slit oriented along the parallactic angle (PA), so their shape can be modified due to differential refraction (Filippenko 1982). Even when the slit is oriented along the PA, contamination due to the light from the host galaxy can produce spurious SEDs. Therefore we have to check the flux calibration of spectra before using them as proxies for the SED. To do the latter, we compute colour indices from the spectra and then we compare them with those obtained using the observed photometry. If the spectrum is well flux-calibrated, then colour differences should be small.
Photometric colour indices at the epoch of the spectra can be computed through the light curve fitting procedure presented in Section 4, while to compute a x1 − x2 colour from a spectrum we use Here, λ is the wavelength in the observer's frame, F λ is the observed SED of the source, S x 1 ,λ and S x 2 ,λ are the transmission functions of the photometric band x1 and x2, respectively, and ZPx 1 −x 2 is the zero point of the colour scale, which can be computed using a star with good spectrophotometric observations. We use the Vega SED published by Bohlin & Gilliland (2004) 4 and the magnitudes published by Fukugita et al. (1996): BVega = 0.03, VVega = 0.03, and IVega = 0.024 mag, and by Cohen et al. (1999): JVega = −0.001, HVega = 0.000, and KVega = −0.001 mag. We adopt the transmission functions given in . Among the SN II spectra available from different sources, we select those: (1) observed in the photospheric phase, and (2) covered by B-and V -band photometry. To check the flux-calibration in the blue and red part of the optical spectra, we compute ∆x 1 −x 2 ≡ (x1−x2) phot −(x1−x2)spec using the B −V and V −I colours, respectively, while for the near-IR spectra we use the J −H and H −K colours, respectively. We also compute the intrinsic B −V colour to represent the shape of the SED. For optical spectra we compute this quantity from dereddened and deredshifted spectra, while for near-IR we compute the intrinsic B −V colour from the photometry (see Section C). values for the collected spectra. For the SN II spectra library, we select spectra with |∆x 1 −x 2 | < 0.1 mag. Finally, we correct spectra for redshift and for Galactic and host galaxy extinction, assuming a Fitzpatrick (1999) extinction curve with RV = 3.1 for both our Galaxy and hosts. Table F5 summarizes the details of the spectra in the library: SN name (Column 1), Galactic colour excess (Column 2), heliocentric redshift (Columns 3), host galaxy colour excess (Column 4), and references for the data (Column 5).

APPENDIX C: C3 LINEARITY
Assuming that a C3 can be well represented by a polynomial fit, the linearity of a C3 can be demonstrated if there is a high fraction of SNe II showing C3s with straight line as optimal polynomial. Due to the scarcity of K-band photometry, we use only BVIJH photometry for this analysis.
With BVIJH photometry set it is possible to define a total of ten colour indices and, therefore, 90 colour-colour plots (i.e., discarding one-to-one correlations). Among them, only 36 combinations give us non-superfluous information, which we analyse for host galaxy colour excess estimation.
Before computing intrinsic C3 slopes, photometry must be corrected for Galactic and host galaxy extinction, and for K-correction. Since we need the prior knowledge of the intrinsic B −V , we need to know in advance the value of host galaxy colour excess. For the latter, we apply zero order correction as prior values. The intrinsic B −V can be computed easily from the relation between the observed and the intrinsic B −V colour, i.e., where (B −V ) obs is the observed B −V . In Sections 5, 6, and 7 we found that R G  For each SN and for each colour combination, we adjust a polynomial fit. The optimum degree is determined using the AIC/BIC and the principle of parsimony. Fig. C1 shows the fraction of SNe that are well represented by a straight line. In 20 of the 36 colour combinations, the number of SNe displaying a linear C3 is over 70 per cent.
Assuming the C3 linearity for different combinations, we compute slopes of all SNe in our set. For each colour combination we compute the median and rms of the C3 slopes. Fig. C2 shows this process for the V −I versus B −V C3s.
Among the spectra available from different sources, we select spectra: (1) of SNe II with explosion epoch constrained within 10 d through photometric information, where for these SNe we assume the midpoint between the last nondetection and the first detection as the explosion epoch (t 0,ln ), (2) spanning a rest-frame wavelength range from 4100 A to 7000Å with a S/N 10Å −1 , and (3) within 40 d since the explosion epoch. We do not include spectra at > 40 d since explosion because in the literature it is less abundant than spectra at earlier epochs (see, for example, Fig. 5 in Gutiérrez et al. 2017) which could bias the age determination toward earlier epochs, and also because at late time the spectral evolution is slower than at early phases, which makes more difficult to accurately determine ages with SNID (Blondin & Tonry 2007). If for a given epoch a SN has several spectra within one day, then we keep that with higher S/N. With the aforementioned constraints, we generate a SNID template library with 242 spectra of 56 SNe II, where each spectrum is corrected by heliocentric redshift. Details of this SN II templates library are summarized in Table F6.  The library has, on average, 6 spectra per day, while almost all the variation is within ±2 rms around the mean.
To create the template library, we run the logwave routine (which is part of the SNID program) with the options w0=3000 w1=8400 nw=1024. This generates SNID spectral templates with a bin in the logarithmic wavelength space of ln(8400/3000)/1024≈0.001, equivalent to 300 km s −1 .
Once the template library is created, the next step is to test how good are the phases since explosion computed with SNID and our new library. To do this, we run SNID with each library spectrum as input, using the avoid option to avoid templates of the same SN. We record all phases and rlap values (which indicate the correlation strength) of the templates found to be similar to the input spectrum and with a redshift within ±0.01. The top panel of Fig. D2 shows the 2D histogram of differences between phases since explosion estimated from last non-detection and from SNID, versus rlap. To correlate rlap with a rms error in phase, we compute the rms of phase differences in bins of width 2 rlap, which is shown at bottom of In general, only one spectrum (e.g., the earliest) is used to estimate explosion epochs with SNID (e.g., Anderson et al. 2014b;Gutiérrez et al. 2017). We expect, however, that including all available spectra of a SN in the analysis will result in a best estimation of the explosion epoch. To explore this possibility, we select all SNe in our library with five or more spectra and perform the following procedure: 1. For each SN we randomly choose one spectrum, computing the explosion epoch (t0,SNID) and t0,SNID − t 0,ln .
2. We compute the median (i.e., the offset) and the rms of the t0,SNID − t 0,ln distribution.
3. We repeat steps 1 and 2 10 2 times, recording the median of the offsets and the rms values. 4. We repeat steps 1-3, but now randomly choosing two and then three spectra per SN as input.  (1989) extinction curve and R V = 1.9 (Pozzo et al. 2006). For that extinction curve, we obtained R h J = 0.402. * (1) IAUC 3532; (2) Thompson (1982); (3) Table D1 shows the result of the aforementioned process, i.e., the offset and the rms as a function of the number of input spectra. Using only one spectrum we obtain a typical rms of 5.0 d, which is similar than the rms of 5.2 d reported by Gutiérrez et al. (2017). We see that using more than one spectrum the rms is reduced down to 4.1 d. The median of the offsets is around -0.6 d, independent of the number of input spectra in the analysis. This offset means that explosion epochs computed with SNID are 0.6 d earlier that those estimated with the non-detection. Hereafter, for the explosion epochs derived with SNID we assume an intrinsic error of 5.0 d when only one spectrum is used, or 4.1 d whether more spectra are available.

APPENDIX E: THE DISTANCE TO NGC 6946
The distance to NGC 6946, host of SN 2004et, was measured with the TRGB method by Tikhonov (2014, hereafter T14), Murphy et al. (2018, hereafter M18), and Anand et al. (2018, hereafter A18), and correspond to µ = 29.39 ± 0.14 mag in the Jang & Lee (2017a) calibration. The PMM Jband distance for SN 2004et obtained in this work (µJ = 28.83 ± 0.12 mag) is in conflict with the TRGB estimation. To investigate the reason for this discrepancy, we compute distances to other two SNe II that exploded in NGC 6946: SN 1980K and SN 2002hh. These SNe have near-IR photometry, but we did not include them into the analysis because they do not have photometry in the five bands we use (i.e., BVIJH).
Using data given in Table D2, we obtain µJ = 28.73 ± 0.18 and 28.77 ± 0.11 mag for SN 1980K and SN 2002hh, respectively, which are consistent with the distance computed with SN 2004et. There is then a tension of at least ∼4 rms between the PMM and the TRGB distance. This discrepancy could be due to: (1) all three SNe II are intrinsically brighter at least 0.56 mag than the SNe II we use for the calibration, or (2) there are issues with the TRGB distances reported by T14, M18, and A18. We noted that the latter two independent studies used almost the same data but obtained significantly different values of the TRGB F814W magnitude (F814WTRGB): 26.00 ± 0.04 mag in M18 and 25.84 ± 0.11 mag in A18. T14 used another set of image data, which is closer to center of the galaxy than those used in M18 and A18, and obtained a lower F814WTRGB value (25.79 ± 0.05 mag). At this moment, the origin of this large discrepancy is unclear. Taking into account this, we safely remove SN 2004et from the calibration and the final sample.

APPENDIX F: TABLES
This paper has been typeset from a T E X/L A T E X file prepared by the author. Table F1. Akaike and Bayesian weights, and evidence ratios for            2456900.5 -1.7 0.0 10 CBET 3964, 29, 24 Column 1: SN names. Column 2: discovery epochs. Column 3 and 4: last non-detection and first detection epochs, respectively, with respect to the discovery epoch. Column 5: spectra phase values with respect to the explosion epoch, which we assume as the midpoint between the last non-detection and the first detection. Adjacent ages are listed in brackets. Column 6: references for data. a Explosion time constraint obtained through polynomial fit to pre-maximum VRI photometry. b C. Feliciano report on the Bright Supernova website (http://www.rochesterastronomy.org/snimages/ ) c (1) Gutiérrez et al. (2017); (2)