Near-IR Type Ia SN distances: host galaxy extinction and mass-step corrections revisited

We present optical and near-infrared (NIR, $YJH$-band) observations of 42 Type Ia supernovae (SNe Ia) discovered by the untargeted intermediate Palomar Transient Factory (iPTF) survey. This new data-set covers a broad range of redshifts and host galaxy stellar masses, compared to previous SN Ia efforts in the NIR. We construct a sample, using also literature data at optical and NIR wavelengths, to examine claimed correlations between the host stellar masses and the Hubble diagram residuals. The SN magnitudes are corrected for host galaxy extinction using either a global total-to-selective extinction ratio, $R_V$=2.0 for all SNe, or a best-fit $R_V$ for each SN individually. Unlike previous studies which were based on a narrower range in host stellar mass, we do not find evidence for a"mass-step", between the color- and stretch-corrected peak $J$ and $H$ magnitudes for galaxies below and above $\log(M_{*}/M_{\odot}) = 10$. However, the mass-step remains significant ($3\sigma$) at optical wavelengths ($g,r,i$) when using a global $R_V$, but vanishes when each SN is corrected using their individual best-fit $R_V$. Our study confirms the benefits of the NIR SN Ia distance estimates, as these are largely exempted from the empirical corrections dominating the systematic uncertainties in the optical.


INTRODUCTION
Since the initial standardization of Type Ia supernova (SN Ia) peak luminosities was employed in the discov-Corresponding author: J. Johansson joeljo@fysik.su.se ery of the accelerated expansion of the Universe (Riess et al. 1998;Perlmutter et al. 1999), estimates of the local value of the Hubble constant from SNe (H 0 Riess et al. 2019) are in tension with the value inferred from the early universe (Planck Collaboration et al. 2020). This tension is a possible sign of new physics or unresolved sources of systematic uncertainty.
Significant work has gone into understanding how to more precisely standardise SNe Ia as distance indicators at optical (visible) wavelengths. The SN Ia optical peak brightness is corrected for lightcurve shape (Phillips 1993) and color (Tripp 1998), and there are now several more elaborated prescriptions for optimising these standardisation procedures (see, e.g., Guy et al. 2007;Burns et al. 2011;Mandel et al. 2011). More recently, additional correction terms aiming at further improving the SN Ia standard candle have also been proposed. One such term accounts for the dependence of the SN Ia luminosity on its host galaxy properties, e.g. stellar mass (Hamuy et al. 1995;Sullivan et al. 2003;Lampeitl et al. 2010;Childress et al. 2013;Betoule et al. 2014;Uddin et al. 2017;Scolnic et al. 2018;Wiseman et al. 2020;Kelsey et al. 2021). These studies all uncover, to various degrees of significance, a "mass step" in the data: after light-curve standardisation, SNe in high-mass galaxies are more luminous than those exploding in low-mass galaxies.
The origin of this mass step is poorly understood, with possible explanations suggesting that it is due to dust in the host galaxies (Brout & Scolnic 2021).
Near-infrared (NIR; 1 < λ < 2.5 µm) observations offer many advantages for standardising SNe Ia (Elias et al. 1985;Meikle 2000). Not only is the NIR less prone to extinction from dust, but SNe Ia are more naturally standard candles at these wavelengths, requiring no or significantly smaller corrections to their peak luminosity to yield similar precision as compared to the optical (Krisciunas et al. 2004;Wood-Vasey et al. 2008;Mandel et al. 2009;Burns et al. 2011;Dhawan et al. 2018a;Burns et al. 2018;Avelino et al. 2019). Theoretical models further corroborate these observations (Kasen 2006;Blondin et al. 2015).
There are already upcoming datasets (e.g. CSP-II, Sweetspot and RAISIN; Phillips et al. 2019;Ponder et al. 2020;Kirshner 2013), and ongoing (e.g. SIRAH, DEHVILS and VEILS; Jha et al. 2019) and future SN Ia programs (e.g. with the Nancy Grace Roman Space Telescope; Hounsell et al. 2018) aiming to take advantage of these properties of SNe Ia and use NIR observations to study dark energy. In this contexts, NIR observations of SNe Ia in the nearby Hubble flow (z > ∼ 0.03) are extremely valuable cosmological tools both as a Hubble flow rung of the local distance ladder and as a low-z "anchor" sample to measure dark energy properties.
However, as Burns et al. (2018) point out, there is a deficit of SNe in low-mass hosts in the current SN Ia NIR data-set and observing an unbiased sample of SNe Ia in the nearby Hubble flow is crucial to test the impact of SN Ia systematics, e.g. extinction from host galaxy dust, on the inferred value of H 0 (Dhawan et al. 2018a;Burns et al. 2018). Moreover, recent works have also claimed evidence for a mass step in the NIR as well (Uddin et al. 2020;Ponder et al. 2020). If indeed present and not accounted for, it will introduce further systematic uncertainties in the NIR SN Ia cosmological analyses.
The main goal of this work is to obtain optical and NIR light curves of an unbiased sample of SNe Ia in the nearby Hubble flow, and together with data from the literature to examine the impact of the host galaxy extinction determination on the claimed correlations between the host stellar masses and the NIR Hubble diagram residuals.
Here we present optical and NIR observations of a new sample of 42 SNe Ia with redshifts out to z ∼ 0.12 and containing 12 SNe in hosts with masses below log (M * /M )=10.
Section 2 presents our sample. Section 3 describes our observations. Section 4 presents our analysis techniques, including spectroscopic classification, light-curve fitting, derivation of the NIR Hubble diagram, and correlations with the host galaxy stellar mass. Section 5 discusses of the results, and Section 6 provides our conclusion.
Throughout this paper we assume flat ΛCDM cosmological model with Ω M = 0.27 and Hubble constant H 0 = 73.2 km s −1 Mpc −1 from Burns et al. (2018).

SUPERNOVA SAMPLE
This work presents 42 new SNe Ia discovered with the intermediate Palomar Transient Factory (iPTF; Rau et al. 2009). We chose targets spanning a wide range of redshifts and host galaxy environments, and acquired optical and NIR follow-up observations for targets with early iPTF detection and classification. These observations are described in more detail in Section 3.
For our analysis, we also include SNe Ia from the literature having both optical and NIR light curves, which we describe briefly here and summarize in Figures 1 and 2. The final photometry of the first stage of the Carnegie Supernova Project (CSP-I) are presented in Krisciunas et al. (2017). Their sample consists of 120 SNe with NIR coverage, z=0.0037 to 0.0835. CfAIR2 (Friedman et al. 2015) is a sample of NIR light curves for 94 SNe Ia obtained with the 1.3m Peters Automated In-fraRed Imaging TELescope (PAIRITEL) between 2005-2011. Barone-Nugent et al. (2012) present J and Hband lightcurves of 12 SNe Ia discovered by PTF in the redshift range 0.03 < z < 0.08. This data was reanalysed by Stanishev et al. (2018), including optical lightcurves. Stanishev et al. (2018) add 16 more SNe with NIR data in the redshift range z=0.037 to 0.183. Furthermore, we include the 6 SNe with UV, optical and NIR lightcurves in Amanullah et al. (2015). Note that some of the supernovae were observed by, e.g., both CSP and CfA (see Friedman et al. 2015, for a comparison), and the total sample size in Figures 1 and 2 refers to the number of unique SNe.

OBSERVATIONS
The follow-up observations were obtained with several different facilities, which are described in the follow-ing sections. For each instrument used, deep reference images were obtained after the supernova emission had faded away. The reference images were subtracted from the science images in order to facilitate the photometry of the SNe, which can otherwise be affected by the light of the host galaxy. Image subtraction was in most cases performed as part of the reduction pipelines, which all utilize implementation of the convolution algorithms presented in Alard & Lupton (1998).

Optical data
During the intermediate Palomar Transient Factory (iPTF) survey, the Palomar 48-inch (P 48) telescope typically delivered g and R-band images. The P 48 image reduction is described by Laher et al. (2014), while the PTF photometric calibration and the photometric system are discussed by Ofek et al. (2012).
Optical follow-up observations were collected using the Palomar 60-inch telescope (P60, BV griz filters), the 2.56 m Nordic Optical Telescope (NOT) and the Las Cumbres Observatory (LCO) in U BV RI and/or griz-bands. The P60 data were reduced using an automated pipeline (Cenko et al. 2006), calibrated against SDSS and the reference images subtracted using FPipe . Similarly, the NOT data were reduced with standard IRAF routines using the QUBA pipeline (Valenti et al. 2011), calibrated to the Landolt system through observations of standard stars and SDSS stars in the field. LCOGT data were reduced using lcogtsnpipe (Valenti et al. 2016) by performing PSF-fitting photometry. Zeropoints for images in the U BV RI filters were calculated from Landolt standard fields (Landolt 1992) taken on the same night by the same telescope. For images in the griz filter set, zeropints were calculated using SDSS magnitudes of stars in the same field as the object.

Near-IR observations
For 37 out of 42 SNe in our sample, we acquired followup observations using the Reionization and Transients InfraRed camera (RATIR). RATIR is a six band simultaneous optical and NIR imager (riZY JH-bands) mounted on the autonomous 1.5 m Harold L. Johnson Telescope at the Observatorio Astronómico Nacional on Sierra San Pedro Mártir in Baja California, Mexico (Butler et al. 2012;Watson et al. 2012;Klein et al. 2012;Fox et al. 2012a).
Typical observations include a series of 80-s exposures in the ri-bands and 60-s exposures in the ZY JH bands, with dithering between exposures. The fixed IR filters of RATIR cover half of their respective detectors, automatically providing off-target IR sky exposures while the target is observed in the neighbouring filter. Master IR sky frames are created from a median stack of offtarget images in each IR filter. No off-target sky frames were obtained on the optical CCDs, but the small galaxy sizes and sufficient dithering allowed for a sky frame to be created from a median stack of selected images in each filter that did not contain either a bright star or extended host galaxy.
Flat-field frames consist of evening sky exposures. Given the lack of a cold shutter in RATIR's design, IR dark frames are not available. Laboratory testing, however, confirms that the dark current is negligible in both IR detectors (Fox et al. 2012b). Bias subtraction and twilight flat division are performed using algorithms written in python, image alignment is conducted by astrometry.net (Lang et al. 2010) and image co-addition is achieved using swarp (Bertin 2010). Figure 3 shows a typical set of images, where blue, green and red shows the field-of-view for i, J and H-band frames, respectively.
For seven SNe Ia in our sample, J and H-band observations were also obtained using other facilities, such as HAWK-I on the 8m Very Large Telescope (VLT) (for iPTF14bbr, 14ddi, 14deb, 14eje and 14fww), VIRCAM on the 4m VISTA telescope (iPTF14fpb) and WIRC on the Palomar 200-inch telescope (iPTF14gnl). These observations were processed with the corresponding instrument reduction pipelines.

Image Subtraction
Image subtraction was performed utilizing the High Order Transform of PSF ANd Template Subtraction (HOTPANTS; Becker 2015). Point sources were selected across the field-of-view (FOV) to calculate the pointspread function (PSF) in each image, either based on classification from SDSS or through manual inspection. Given the relative paucity of bright point sources in most fields (particularly in the NIR), the PSF was held fixed across the FOV.
The calculated PSFs were utilized to perform PSFmatched photometry on the resulting subtracted images, yielding measurements of the instrumental magnitude of the supernova in each epoch: Uncertainties and upper limits were determined by inserting false sources of varying brightness into the RATIR images and repeating the identical process of image subtraction and PSF-matched photometry.

Photometric Calibration
Photometric calibration of the RATIR data was performed following the process outlined in Ofek et al. (2012). To calculate color and illumination terms, we selected fields with coverage from both SDSS (optical) and UKIDSS (NIR), and obtained photometry for stars (i.e., objects classified as point sources in SDSS) with r-band magnitudes between 14 and 18 (with additional flagging for saturation). We measured instrumental magnitudes via PSF-matched photometry for these calibration stars as above.
As a first pass, we calculate a zero-point for each image with no additional corrections (e.g., color and illumination terms). We removed nights with large scatter in the zeropoint (RMS ≥ 0.10 mag) or individual stars that were clear outliers in the fits (determined via visual inspection).
We then performed a least squares fit using the remaining nights/stars to the following equation for each filter f : where ZP is the zero-point, CT f,i,j is the color term, m i and m j are the filters used for the color correction, and C illum is an illumination correction term accounting for PSF variations depending on the position on the detector.
The color and illumination terms were held fixed for all observations in a given filter, while the zero-point term was allowed to vary freely in each image. The resulting best-fit color terms and zero-point RMS are shown in Table 1. The zero-point RMS is typically ∼0.03 mag for the RATIR r to H-band.
For fields with SDSS and UKIDSS coverage, calibrated supernova magnitudes were calculated using Equation 2. For fields lacking SDSS coverage, we used photometry from Pan-STARRS1 Data Release 2 (Magnier et al. 2020), which is in a photometric system close to SDSS (Tonry et al. 2012). For fields lacking UKIDSS coverage we used 2MASS (Skrutskie et al. 2006) and the transformation from Hodgkin et al. (2009) to calibrate the Y -band RATIR data (Y = J +0.50×(J −H)+0.08).

ANALYSIS
Some of the SNe presented here have previously been published in separate papers: • Optical and NIR light curves and spectra of iPTF13abc (SN 2013bh) were presented and analysed in Silverman et al. (2013). It is a near identical twin to the peculiar Ia SN 2000cx.
• UV, optical and NIR light curves and spectra of iPTF13asv (SN 2013cv) were presented in Cao   Weyant et al. (2018). iPTF13asv shows low expansion velocities and persistent carbon absorption features after the maximum, both of which are commonly seen in super-Chandrasekhar events, although its light curve shape and sharp secondary near-IR peak resemble characteristic features of normal SNe Ia.
• Optical light curves and high-resolution spectra of iPTF13dge were presented in Ferretti et al. (2016), and NIR light curves in Weyant et al. (2018). The light curves are compatible with a normal SN Ia with little reddening, and no definite timevariability could be detected in any absorption feature of iPTF13dge.
• UV, optical and NIR observations of iPTF13ebh from the CSP-II collaboration were presented in Hsiao et al. (2015). iPTF13ebh can be categorized as a "transitional" event, on the fast-declining end of normal SNe Ia, showing NIR spectroscopic properties that are distinct from both the normal and subluminous/91bg-like classes.
• iPTF14atg is a subluminous peculiar SN similar to SN 2002es. It displayed strong, declining ultraviolet emission shortly after explosion. Spectra together with UV, optical and NIR photometry have been extensively analysed in Cao et al. (2015); Kromer et al. (2016).
• UV and optical photometry and spectra of the 1999aa-like SN iPTF14bdn were presented in Smitka et al. (2015). The rapid, near-linear rise, the non-evolving blue colors, and strong absorption from ionized carbon, are interpreted to be the result of either vigorous mixing of radioactive-Ni in the SN ejecta, or ejecta interaction with diffuse material, or a combination of the two.
• iPTF17lf was reddened, spectroscopically normal SN Ia, discovered during a wide-area (2000 deg 2 ) g and I-band survey for "cool transients" as part of a two month extension of iPTF (Adams et al. 2018).
For the other SNe included in our sample (except iPTF14ale, which has no spectroscopic classification), we run the SuperNova IDentification code (SNID Blondin & Tonry 2007) on the spectra (to be presented in a separate paper). For SNe iPTF13s, iPTF13ddg, iPTF13efe, iPTF14bpz, iPTF14fpb we rely on redshift estimates based on the SN spectral features using SNID. Furthermore, for iPTF13anh, iPTF13asv, iPTF13azs, iPTF13crp and iPTF13dkx, we determine the redshifts from narrow host galaxy lines in the SN spectra. For our single spectrum of iPTF14apg, observed 5 days before peak brightness, SNID gives a best match to SN 2004dt at z = 0.088 ± 0.004, consistent with the spectroscopic redshift of the nearest galaxy. 4.1. Host galaxies Figure 12 shows cut-out images from the SDSS and PanStarrs surveys, centered on the SN positions. Most SNe can easily be associated with their hosts, while some cases are ambiguous, including: • iPTF14apg: nearest galaxy is SDSS J123758.69+082301.5 with a spectroscopic redshift z=0.08717, separated by 51 , corresponding to a projected distance of 79.4 kpc.
For the literature sample, we note that SNe PTF10hmv, PTF10nlg and PTF10qyx from Barone-Nugent et al. (2012) have ambiguous hosts.
We estimate the host galaxy stellar mass, M * , using the relationship published in Taylor et al. (2011), We use g and i-band magnitudes from SDSS (or PanStarrs when no SDSS photometry was available), corrected for the Milky Way (MW) extinction. M i is the absolute magnitude in the i-band. Table 2 lists the redshifts and coordinates of the SNe in our sample, together with their likely host galaxies and our estimates of the host galaxy stellar mass. Our mass estimates are consistent with those of (Neill et al. 2009;Chang et al. 2015;Burns et al. 2018), but systematically higher by ∼0.2-0.3 dex than the estimates from Ponder et al. (2020) and Uddin et al. (2020), who employ a more sophisticated SED fitting (Fig. 4). However, for consistency when comparing stellar masses between our sample and the CSP, CfA and literature sample, we choose to use our estimates for the combined analysis.

Light curve and host galaxy extinction fitting
We use the SNooPy light curve fitting package developed for the CSP sample (Burns et al. 2011(Burns et al. , 2014(Burns et al. , 2018 to analyze the light curves the SNe in our sample, including the light curve of the literature sample. To find the time of maximum, T max , and color-stretch parameter, s BV , and the observed rest-frame peak magnitudes 1 of the SNe, the SNooPy max model was fitted to the light curves. An example fit is shown in the upper panel of Figure 5 and the derived light curve parameters are given in Tables 3 and 4. To derive the host galaxy extinction we use the more elaborated color model.
Upper panel shows the max model fit to iPTF13dkx, a representative SN from our sample at "moderate" redshift and NIR coverage around peak brightness. Lower panel shows the inferred color excess (normalized with respect to V −band) and the best-fit extinction parameters. the host galaxy extinction taking into account the dependence of SN Ia intrinsic colors on s BV (Burns et al. 2014). It uses parametrized dust extinction laws to calculate the total-to-selective extinction ratio R X in any filter X, as a function of R V and E(B − V ) host by the means of synthetic photometry (see Burns et al. 2011). As R V controls the wavelength dependence of the extinction and the host-galaxy color excess E(B − V ) host the amount of the extinction, with observations over a broad range of frequencies it is in principle possible to fit independently for R V and E(B − V ) host , which are otherwise correlated. In our analysis we used Cardelli et al. (1989) and O'Donnell (1994) extinction law. For full details on the color model the reader is referred to Burns et al. (2014Burns et al. ( , 2018. We performed two fits for the extinction. First, R V = 2.0 was assumed for all SN hosts and only E(B − V ) host was fitted. The value R V = 2 corresponds to our sample average (weighted average R V = 1.9, σ R V = 0.8) and is close to values commonly found in many SN Ia cosmological analyses, which commonly employ single R V . Second, both E(B − V ) host and R V were fitted. This is possible because SNe Ia show a small intrinsic color dispersion across optical to NIR bands and the wavelength leverage provided by including NIR observations. Nevertheless, when E(B − V ) host approaches zero (or rather the level of scatter in the intrinsic color, σ E(B−V ) ∼ 0.06 mag), the leverage to get meaningful constraints on R V decreases. The results from the second fit are shown in Table 4. Figure 5 lower panel shows an example of the inferred color excess and the best-fit extinction parameters.
It is a long-standing issue that SN analyses have yielded "unusually" 2 low R V values. This is seen both when minimizing the Hubble residuals using a global R V for cosmological samples and for detailed studies of individual, highly-extinguished SNe (e.g. SNe 2006X and 2014J; Burns et al. 2014;Amanullah et al. 2014). We stress that we only use the observed colors to constrain E(B − V ) host and R V , since determining extinction by minimizing Hubble residuals can lead to a bias (Burns et al. 2018;Uddin et al. 2020).

NIR Hubble diagram
To construct the Hubble diagrams, the distance modulus for filter X, µ X , was computed as: where P N X (s BV − 1) is the 2-nd order polynomial luminosity-decline-rate relation from Burns et al. (2018) and R X,BV is the total-to-selective absorption coefficient for filter X computed from R V and E(B − V ) host . Here, we impose that R V > 0 and do not correct for dust extinction objects with E(B − V ) host < 0, i.e. intrinsically blue objects. Figure 6 shows the resulting J and H-band Hubble diagrams for our optical+NIR SNe Ia compilation, including 40 of SNe from our sample presented. iPTF14apg and iPTF14atg, are not included here, as we do not include spectroscopically peculiar SNe Ia (03fg, 06bt, 02eslike nor Iax SNe) in the analysis. Furthermore we apply a set of cuts on the redshift, stretch and color excess distribution on our sample, such that we include only SNe with z CM B > 0.01, s BV > 0.5, E(B − V ) host < 0.5 mag, and E(B − V ) MW < 0.2 mag (corresponding to typical sample cuts used in other cosmological analyses, e.g. using SALT2 parameters −0.3 < c < 0.3 and −3 < x 1 < 3). The solid lines show the best-fit Hubble lines and the dashed lines indicate the scatter expected due to peculiar velocities v pec = 300 km/s. The RMS scatter in the Hubble residuals for the combined sample, after the cuts, is σ HR,J =0.19 mag (165 SNe) and σ HR,H =0.21 mag (152 SNe), for J and H respectively. The scatter in J and H does not decrease significantly when using individual best-fit R V instead of a global R V .
We note an offset of 0.20 ± 0.05 mag when comparing the Y -band peak magnitudes to the CSP-I sample (also seen when comparing individual Y -band light curves of SNe observed simultaneously by RATIR and CSP-II, private communication). We thus add 0.20 mag to the Y -band magnitudes listed in Table 3 for the Hubble diagram analysis.

Correlations with host galaxy stellar mass
Having SN host galaxy stellar masses determined in Sect. 4.1 and color and stretch corrected distances from Sect. 4.3, we can begin to look for correlations.
In Figure 7 we show how our derived color stretch and color excess correlate with host stellar mass. Similar to conclusions reached in previous studies (Sullivan et al. 2011;Childress et al. 2013), we find that low-mass galaxies tend to host SNe with higher stretch (s BV > 0.8) with moderate extinction (E(B − V ) host < ∼ 0.25 mag), while high-mass galaxies also host highly reddened SNe and fast-declining SNe.
Following Stanishev et al. (2018) and references therein, we fit the probability density function (PDF) of the computed color excesses for the entire sample, using an exponentially modified Gaussian distribution with a mean c 0 and standard deviation σ c and exponent relaxation parameter τ . We find values c 0 =0.02 mag and σ c =0.06 mag and τ =0.14. We interpret the Gaussian component as a residual scatter due to intrinsic color variations.
Previous analyses (e.g. Sullivan et al. 2011;Betoule et al. 2014) typically split the sample at M split = 10 10 M , which seems to be an "astrophysically reasonable" choice given the fairly distinct difference between the stretch and color excess distributions below and above log ( spective sample or based on some information criterion that maximizes the likelihood (Uddin et al. 2020;Ponder et al. 2020;Thorp et al. 2021). We choose to split our sample at M split = 10 10 M , as our fiducial case. Despite adding more SNe in low-mass galaxies from our sub-sample and e.g. the Barone-Nugent et al. (2012) sub-sample, the distribution of host stellar mass for our sample is still skewed towards higher log (M * /M ). For the combined sample, the median log (M * /M )=10.50.
If we look at the observed distribution of best-fit R V values (Fig. 8), we find a weighted average R V = 2.2 (σ R V = 0.9) for log (M * /M )< 10.0 and R V = 1.7 (σ R V = 0.8) for log (M * /M )> 10.0 host galaxies. The weighted average value of R V for the whole sample is R V = 1.9 (σ R V = 0.9). Here, we are not including R V estimates for SNe with color excesses close to the level of intrinsic color scatter E(B − V ) host < σ c ∼ 0.06 mag (where we typically find artificially low R V , albeit with large error-bars) nor for highly extinguished SNe E(B − V ) host > 0.5 mag (which are well fit by R V values ranging from 1.1 to 2.7, but the distribution is likely to be observationally biased towards finding SNe with low R V ).
In order to test the hypothesis that the distributions of R V in the low and high stellar mass bins are statistically compatible from being drawn from the same underlying distribution we perform a Kolmogorov-Smirnov (K-S) test. The K-S test yields a p-value of only 0.015, hence suggesting that the distributions are significantly different, with more than 95% confidence level.
Even though the weighted mean values are statistically consistent with the global mean R V value, we stress that the wide, non-gaussian, probability distributions are different.
For the low-and high-mass bins, we compute a weighted mean and standard deviation of the Hubble residuals for BV griY JH filters (horizontal black lines in Fig. 9 and 11). For each filter X, we will refer to any Hubble residual offset between the two bins as a "massstep", ∆ HR (X). To further investigate the behaviour of the Hubble residuals, we divide the sample into five mass bins (orange symbols in Fig. 9), to see if there are any additional effects towards the edges of the host mass distribution. Following Uddin et al. (2020), we also fit a slope to the Hubble residuals as a function of host mass (yellow lines in Fig. 9 and 11) using Orthogonal Distance Regression (ODR).
We find that using a global value R V = 2.0 (close to the average R V for the entire sample) for all SNe in low-and high-mass host galaxies, we reproduce a significant (∼ 2σ) mass-step in optical BV gri-band Hubble residuals ∆ HR ∼ −0.07 ± 0.03 mag, while for NIR JH-bands there is no significant mass step (∆ HR (J) = −0.021 ± 0.033 mag and ∆ HR (H) = 0.020 ± 0.036 mag), shown as red symbols in Fig. 10 (left panel). A similar trend is seen when fitting a slope to the Hubble residuals as a function of host mass. For optical BV gri-bands we find a slopes of ∼ 0.06 ± 0.02 mag/dex, while in the NIR the slopes are smaller (−0.027±0.016 mag/dex and −0.005±0.018 mag/dex in J and H, respectively) shown as red symbols in right panel of Fig. 10.
When correcting each SN individually by their bestfit R V × E(B − V ) host we see no significant mass-step or mass-slope in the Hubble residuals, across the optical and NIR bands (blue symbols in Fig. 10).
This result seems valid when changing the cuts on z, s BV and E(B − V ) host , and perhaps more importantly the choice of M split . Choosing M split = 10 10.5 M (the median stellar host galaxy mass of our SNe compilation), we do see a (non-significant) mass-step across optical and NIR ∆ HR ∼ −0.04 ± 0.03 mag. We find that our results are in line with with Brout & Scolnic (2021), who modelled host galaxy reddening as separate Gaussian distributions for galaxies below and above log (M * /M )=10. They found that SNe in lowmass hosts, the average R V = 2.75 ± 0.35, whereas for SNe in high-mass hosts, R V = 1.5 ± 0.25, with both sub-samples having a wide distributions σ R V = 1.3. This is in fair agreement with Salim et al. (2018), who find that on average, dusty, high-mass quiescent galaxies have lower R V values ( R V = 2.61), whereas low-mass star forming galaxies tend to have higher values for R V ( R V = 3.15). Uddin et al. (2020) found nominal evidence for a consistent mass-step in both the optical and NIR using the CSP-I sample (∆ HR,J = −0.103 ± 0.050 mag, and ∆ HR,H = −0.097 ± 0.047 mag) using similar cuts on the sample, although including SNe with z < 0.01 and having the mass-step located at log (M * /M )=10.5 (shown as gray dashed lines in Fig. 10). We can not fully reproduce the NIR mass-step reported by Uddin et al. (2020), even if we use their host masses and the same best-fit extinction. We note that they do use updated Phillips-relations, correcting for stretch using more flexible spline functions calibrated using unpublished data from CSP-II (C. Burns and S. Uddin, private communication), but we do not expect this to be result in the observed differences in our plots. Ponder et al. (2020) also report a H-band massstep ∆ HR,H = −0.18 ± 0.05 mag (mass-step located at log (M * /M )=10.44) using a compilation of 99 SNe from the literature. However, after removing two outliers, the step reduces to ∆ HR,H = −0.10 ± 0.04 mag. It is unclear if their results are due to the lack of NIR stretch corrections and/or color-corrections.
Recently, Thorp et al. (2021) analyzed optical (griz) lightcurves of 157 nearby SNe Ia (0.015 < z < 0.08) from the Foundation DR1 dataset using the BayesSN lightcurve fitter (Mandel et al. 2011(Mandel et al. , 2020. When split-  Table 2) for the same cuts on E(B − V ) host < 0.5 mag and sBV > 0.5, but for log M split = 10.5 and including SNe at z < 0.01. The dotted line shows the mass-step in Thorp et al. (2021). ting their sample at log (M * /M )=10, they find R V = 2.84 ± 0.31 in low-mass hosts, and R V = 2.58 ± 0.23 in high-mass hosts. They conclude that these values are consistent with the global value of R V = 2.61 ± 0.21, estimated for the full sample, and can not be an explanation of the mass step. After corrections, their resulting mass-step is ∆ HR = 0.050 ± 0.022 mag (shown as gray dotted line in left panel of Fig. 10).

DISCUSSION
The relation between SN Ia luminosity and host galaxy properties is of great interest for SN Ia progenitor studies as well as for cosmology, as a third empirical correction (Sullivan et al. 2010;Lampeitl et al. 2010;Kelly et al. 2010). While the correlation has been seen nearly ubiquitously across different samples (see e.g. Brout et al. 2019;Smith et al. 2020), there is a significant debate about the physical origin of this relation. There has been speculation that the mass-step may be driven by the age of the stellar population, metallicity or star formation rate (Sullivan et al. 2010;Gupta et al. 2011;D'Andrea et al. 2011;Hayden et al. 2013;Childress et al. 2014;Rigault et al. 2020;Rose et al. 2020). However, it is possible that this correlation arises due to different SN Ia environments in different host galaxy types.
While most studies pertaining to SN luminosity and host galaxy correlations in the optical, recently there have been reports of the possible detection of the massstep in the NIR wavebands (Ponder et al. 2020;Uddin et al. 2020). In this study, we exploit the multiwavelength, well-sampled lightcurves of the SNe in our sample and compute the mass-step/slope after stretch and color corrections, fitting the R V parameter individually for each SN. We find that when fitting the R V value for each SN, we see a mass step consistent with zero in all filters from B to H−band (see Figure 10). However, when fixing the R V value to the sample average, as is done in previous studies, we find that there is a mass step of ∼ 0.07 -0.1 mag in the optical (BV gri) while no significant step is seen in the NIR (Y JH).
Since the free R V case yields mass steps that are consistent with zero, it is likely that the origin of the mass step is due to variations of dust properties in the interstellar medium of the host galaxies. However, more detailed studies would be required to rule out intrinsic effects. In fact, there are indications that there are two SN populations, having different SN ejecta velocities and different intrinsic colors, which also trace the host galaxy stellar mass (see e.g. Polin et al. 2019;Siebert et al. 2020;Pan 2020). Childress et al. (2013) and Gonzalez-Gaitan et al. (2020) simulate the effect of having separate colorluminosity corrections for low-and high-mass galaxies. They find that multiple free color-luminosity slope parameters may explain away the mass-step, suggesting that the origin of mass-step is a difference in intrinsic color-luminosity relation (β int c int ) of two SN populations found in galaxies with different masses as opposed to different dust properties.

CONCLUSION
Many studies in the literature have suggested the need for an additional standardization parameter for SNe Ia, beyond the lightcurve shape and color. In particular, firm evidence has been put forward for a correlation between residuals in the Hubble diagram at optical wave-lengths, and the host galaxy stellar mass (Sullivan et al. 2010;Lampeitl et al. 2010;Kelly et al. 2010;Gupta et al. 2011;D'Andrea et al. 2011;Hayden et al. 2013;Childress et al. 2014;Rose et al. 2020). The underlying cause of these correlations is not completely understood, with some suggestions that this could be due to correlation with age/metallicity of the underlying stellar population, however, there is also evidence pointing to this correlation arising from differences in dust properties of the SN hosts.
Our work differs from previous studies in that we use a sample of SNe Ia, found in the untargeted PTF/iPTF surveys, which adds a significant number of low-mass host galaxies. Furthermore, our data-set includes multiwavelength follow-up observations, including near-IR, which allows us to infer the total-to-selective extinction parameter, R V , for each SN individually. This is motivated by the findings in e.g. Amanullah et al. (2015) and Burns et al. (2018), suggesting that the wavelength dependence of dimming by host galaxy mass varies between SNe, making the use of a single value of R V questionable. Using a parameterized extinction relation by CCM, we fit for both the color-excess, E(B−V ), and R V using SNooPy color model fits of the multi-band data. In other words, our estimate of the extinction along the line-of-sight of each SN is not derived from the Hubble residuals. When examining the Hubble residuals, we do not find a significant correlation with stellar mass at optical or NIR wavelengths.
If we, instead, assume a single value for R V to correct all SNe fitting only the color excess, we recover the "mass-step" in optical filters. In the NIR, we find no significant dependence on the stellar mass, independently of how R V is measured, i.e., individually or globally. This is consistent with the interpretation made by (Brout & Scolnic 2021), that the mass-step is likely caused by differences in dust properties of the low-and high-mass SN host galaxies.