Photometric study of the late-time near-infrared plateau in Type Ia supernovae

We present an in-depth study of the late-time near-infrared plateau in Type Ia supernovae (SNe Ia), which occurs between 70-500 d. We double the existing sample of SNe Ia observed during the late-time near-infrared plateau with new observations taken with the Hubble Space Telescope, Gemini, New Technology Telescope, the 3.5m Calar Alto Telescope, and the Nordic Optical Telescope. Our sample consists of 24 nearby SNe Ia at redshift<0.025. We are able to confirm that no plateau exists in the Ks band for most normal SNe Ia. SNe Ia with broader optical light curves at peak tend to have a higher average brightness on the plateau in J and H, most likely due to a shallower decline in the preceding 100 d. SNe Ia that are more luminous at peak also show a steeper decline during the plateau phase in H. We compare our data to state-of-the-art radiative transfer models of nebular SNe Ia in the near-infrared. We find good agreement with the sub-Mch model that has reduced non-thermal ionisation rates, but no physical justification for reducing these rates has yet been proposed. An analysis of the spectral evolution during the plateau demonstrates that the ratio of [Fe II] to [Fe III] contribution in a near-infrared filter determines the light curve evolution in said filter. We find that overluminous SNe decline slower during the plateau than expected from the trend seen for normal SNe Ia


INTRODUCTION
Although Type Ia supernovae (SNe Ia) are widely used as cosmic distance indicators (Riess et al. 1998, Perlmutter et al. 1999), there is still debate about their explosion mechanisms and the nature of their progenitors (see Hillebrandt et al. 2013, Maoz et al. 2014, Ruiter 2020, Jha et al. 2019 for comprehensive reviews).It is generally accepted that SNe Ia originate from the thermonuclear explosions of carbon-oxygen white dwarfs (CO WDs).The CO material burns to iron-group elements, and the radioactive decay of 56 Ni → 56 Co ( 1/2 = 6 d) powers the early light curves of SNe Ia.Around 60 d post explosion the dominating radioactive decay chain shifts to the decay from 56 Co → 56 Fe, with a longer half life of 78 d.
In the nebular phase (phase > 150 d), the outer layers of the SN ejecta become transparent and the inner regions of the ejecta become visible.Late-time spectroscopy can be used to search for hydrogen, which would point towards a single-degenerate scenario (Hamuy et al. 2003, Mattila et al. 2005, Leonard 2007, Shappee et al. 2013, Silverman et al. 2013, Lundqvist et al. 2013, 2015, Maguire et al. 2016, Kollmeier et al. 2019, Graham et al. 2019, Prieto et al. 2020).Late-time spectra can also be used to constrain the amount of stable nickel vs. unstable material, which can be compared to predictions from explosion models where a larger ratio points to burning at higher central densities, implying a larger progenitor mass (Mazzali et al. 2015, Botyánszki & Kasen 2017, Maguire et al. 2018, Flörs et al. 2020).
The majority of studies at all epochs focus on the optical, because SNe Ia are brightest at these wavelengths and there are many optical instruments available.Studying the near-infrared (NIR,  > 0.8 µm) is more difficult, but is beneficial because SNe Ia are better standard candles in this wavelength range (Elias et al. 1981, 1985, Krisciunas et al. 2004, Wood-Vasey et al. 2008, Barone-Nugent et al. 2012, Johansson et al. 2021, Jones et al. 2022, Galbany et al. 2022, Müller-Bravo et al. 2022).They are less impacted by extinction, with inferred distance estimate root mean square values dropping by 2-4 even before any corrections are applied to the light curves (Avelino et al. 2019).
While most NIR studies focus on light curve standardisation around peak (−5-40 d), observations around 70-600 d are very important for understanding the evolution of the ejecta.Our understanding of the NIR evolution of SNe Ia at late times has evolved significantly during the past four decades.Axelrod (1980) predicted an "IR-catastrophe" -a shift from optical and NIR emission lines to fine-structure iron lines in the mid and far-IR that occurs at around 450 d due to an onset of thermal instability that causes a dramatic temperature change from ∼3000 K to ∼300 K.However, the resulting sharp decline in the optical and NIR light curves has never been observed in SNe Ia.Graur et al. (2020) found instead that SNe Ia reach a plateau in the  and  bands starting at ∼150 d and lasting for approximately a year.The presence of a plateau in the NIR was first predicted by Fransson et al. (1996), where it was linked to the onset of the IR-catastrophe.Fransson et al. (1996) suggested the flattening of the J band is due to the shift from emission of [Fe ] at ∼ 5000 Å to emission of [Fe ] at 1.257 µm and 1.644 µm, which is supported by the NIR evolution of SN 2014J presented by Diamond et al. (2018).Updated spectral models by Fransson & Jerkstrand (2015) showed that a redistribution of ultra-violet (UV) emissivity increases the flux in the optical and NIR, circumventing the "IR-catastrophe".
This flux-redistribution behaviour is reminiscent of the rebrightening seen in the NIR around 30 d past maximum, also called the secondary maximum.This feature is caused by sharp peaks in the emissivity of iron/cobalt gas at certain temperatures, which are near an ionisation edge (Kasen 2006).The dependence of the emissivity on temperature explains why for subluminous and cooler 1991bglike SNe Ia (Filippenko et al. 1992, Leibundgut et al. 1993, Turatto et al. 1996, Taubenberger 2017), the secondary maximum is shifted to earlier phases, often causing it to blend with the primary maximum.The strongest peak in the NIR emissivity occurs at ∼7000 K and represents the ionisation edge between doubly ionised to singly ionised iron, when the ejecta becomes very efficient at redistributing flux from the UV to longer wavelengths, which leads to the re-brightening in the NIR during the secondary maximum.Another peak exists at the ionisation edge between singly ionised and neutral iron at ∼2500 K, which may coincide with the onset of the NIR plateau at ∼150 d.Diamond et al. (2018) presented NIR spectra of SN 2014J during the plateau phase, which demonstrated a decrease in the strength of [Fe ] features in favour of [Fe ] features, but no [Fe ] features.Sollerman et al. (2004) and Graur et al. (2020) agreed that the scattering of UV photons to longer wavelengths is the most likely cause of the NIR plateau.
The end of the plateau at ∼500 d is not yet understood, although a tentative detection of [Fe ] by Graur et al. (2020) suggested a third shift in the dominant ionisation state of iron.Tucker et al. (2022) also identified strengthening features in the optical after the end of the plateau that could be attributed to [Fe ].Graur et al. (2020) tentatively suggested that the plateau does not extend to the   -band, and that the plateau in the H band is comprised of two distinct branches.No theoretical explanation for this bimodal behaviour was offered, although a correlation between the peak magnitude and the magnitude of the plateau would align well with the idea that the plateau is caused by a similar mechanism as the secondary maximum.
In this paper we extend the sample of SNe Ia with NIR photometry on the late-time plateau to 24 SNe Ia, and use the additional data to confirm the absence of the plateau in the K  band, and test whether the magnitudes in the H-band plateau consist of two distinct branches as suggested by Graur et al. (2020).In Section 2, we introduce the sample of nearby SNe Ia and describe the late-time NIR photometry and spectroscopy included in this paper.Fitting methods implemented in this work, as well as the radiative transfer models of SNe Ia in the nebular phase that we compare to our data are described in Section 3. We present our results in Section 4 and discuss their implications in a theoretical context in Section 5. Finally, we summarise and conclude in section 6.

DATA
In Section 2.1, we present the sample of nearby SNe Ia used in this paper.We describe the spectra and photometry included in this paper, which is a combination of data pulled from the literature and new data, in Section 2.2.

Sample of nearby SNe Ia
Our sample consists of new data, as well as data that have previously been published, totalling 24 SNe Ia.Of these, 20 are classified as normal SNe Ia, two are classified as 91T-like (SNe 2000cx and 2021wuf), and two are classified as transitional objects (SNe 2004eo and2012ht).The NIR photometry of six SNe Ia (SNe 2020ees, 2020uxz, 2021jad, 2021pit, 2021wuf, and 2021aefx), and XShooter spectra of two SNe Ia (SNe 2016hvl and 2017cbv) are presented for the first time in this paper.XShooter spectra of three SNe (SNe 2012cg, 2012ht, and 2013aa) published by Maguire et al. (2013), and spectra of four SNe (SNe 2012fr, 2013cs, 2013ct, and 2013dy) published by Maguire et al. (2016) are also included.We performed synthetic photometry on these spectra to extract NIR photometry as described in Section 2.2.2.We include NIR photometry for ten SNe Ia presented by Graur et al. (2020) (SNe 2000cx, 2001el, 2004eo, 2011fe, 2012ht, 2013dy, 2014J, 2017erp, 2018gv, and 2019np, which were originally published by Krisciunas et al. 2003, Sollerman et al. 2004, Pastorello et al. 2007, Stritzinger & Sollerman 2007, Sand et al. 2016, Shappee et al. 2017, Burns et al. 2018).We note that the NIRI observations of SN 2020uxz were taken in K rather than K  .The early observations of SN 2001el (Krisciunas et al. 2003) are taken in a mixture of K and K  , but here we only include the data taken in K  .The late-time observations from Stritzinger & Sollerman (2007) are taken exclusively in K  .We never mix K and K  data when performing fits, as will be discussed further in Section 4.An overview of the sample is presented in Table A1.
Observing the NIR plateau is difficult because SNe Ia are inherently fainter in the NIR compared to the optical and by 150 d they have faded by ∼6 magnitudes relative to peak.Consequently, all the SNe Ia in our sample are nearby, have  B max < 15 mag (with the exception of SN 2016hvl which has  B max =15.4 mag), and are offset from their host galaxies to reduce host contamination.We note that the last criterion is one potential source of bias in our sample (Wang et al. 2013).
All distance moduli and uncertainties were taken either from the literature where available, or calculated from redshift-independent distances provided by the NASA/IPAC Extragalactic Database (NED). 1 The data sources for the distance moduli are summarised in Table A1.
All the photometry measurements were obtained using the package (Brennan & Fraser 2022). 3The data were calibrated using the 2MASS catalog in the Vega magnitude system.Since none of the sources are in very crowded fields, we implemented aperture photometry for the whole sample.As an additional test, point-spread function (PSF) photometry was performed where possible and compared to the aperture photometry.The aperture and PSF magnitudes were consistent within the uncertainties for all measurements.All sources are bright and far removed from their host galaxy and therefore template image subtraction was not required.Background surface fitting failed for a subset of the sample due to the noisy nature of the NIR images, so we reverted to local background fitting for the whole sample.
No S-corrections were applied to the HST photometry because no synchronous J/F125W or H/F160W data were available.We estimate the systematic offset between the filters by performing synthetic photometry on all the XShooter spectra in our sample for the J, H, F125W and F160W bands.On average, we find that the F125W photometry is 0.3 mag fainter than the J band, and F160W is 0.4 mag fainter than the H band.We do not correct for these offsets but any HST photometry is highlighted in Fig. 1 and the reader should note that these points are expected to be fainter than the corresponding ground-based filters.
Finally, all the data were corrected for Milky Way extinction using the dust map provided by Schlafly & Finkbeiner (2011) and the Python module (Green 2018).The photometry was not corrected for host galaxy extinction because all the SNe are well separated from their host, and NIR photometry is minimally impacted by extinction.All the new NIR light curves, together with the NIR light curve data presented by Graur et al. (2020), are shown in Fig. 1.

Spectroscopy
We include 12 mid-resolution spectra of eight SNe Ia obtained using XShooter.Eight of these spectra were previously presented by Maguire et al. (2013Maguire et al. ( , 2016)).Four spectra are published here for the first time and were reduced using the same method described by Maguire et al. (2016).Due to the relatively high spectral resolution of XShooter (∼35 km s −1 ), host galaxy features were easily identified and removed in the reduction process (Maguire et al. 2016).We do not expect to see any contribution from a potential companion star since the remnant models presented in (Pan et al. 2012) predict that the contribution will be very faint relative to the SN at these phases.The spectral response of XShooter is relatively stable with a relative flux uncertainty across the spectrum of 5 per cent (Vernet et al. 2011).The three arms of the XShooter spectrograph were firstly combined using their overlap wavelength regions with small scalings in their flux levels.Next, the spectra were flux calibrated using photometry from Las Cumbres Observatory Global Telescope Network (LCO; Brown et al. 2013) if possible, or alternatively, using photometry performed on stars in the Sloan Digital Sky Survey (SDSS) or Pan-STARRS1 (PS1) field of the acquisition image.As a last resort, zeropoints were taken from XShooter.
The XShooter acquisition image of each SN was used to estimate the magnitude of the SN and comparison stars in the field of the SN.The spectra of SNe 2012cg, 2012ht, and 2013ct were calibrated by comparison of the companion stars to catalogue magnitudes from the SDSS Data Release 10 ( Ahn et al. 2014), and the spectrum of SN 2016hvl was calibrated by comparison to the PS1 Data Release 2 catalog (Flewelling et al. 2020).The spectrum of SN 2017erp was calibrated to LCO photometry of the SN itself taken at similar phases to the spectral observation.The spectra of SNe 2012fr, 2013aa, and 2013cs were calibrated using the XShooter zero-points because no coeval SN photometry nor catalogue magnitudes from SDSS or PS1 were available.In these cases, the tabulated XShooter zero-point was used, resulting in a larger uncertainty.The uncertainty was estimated by comparing the magnitudes obtained using the zero-point method for SNe that also had catalogue measurements.These were found to be <0.5 mag, which we set as the conservative uncertainty of the magnitudes estimated using the zero-point method.The different sources for flux calibration result in a large range of uncertainties.A summary of the flux calibrations is presented in Table 3.
XShooter spectra extend from 5000 Å to 25000 Å, but in some cases the spectrum is very noisy at the far red end.We excluded spectra with spurious flux values at the red edge of the detector by visual inspection.We used (Barbary et al. 2022) to integrate across the J, H, and K  2MASS bandpasses to obtain synthetic photometry (see Table 3 and Fig. 1).

METHODS
In Section 3.1, we describe how we fit the NIR data to determine if there is a plateau, and how we derive the average magnitude and decline rate of the plateau.In Section 3.2, we describe SALT3 fits performed on the optical light curves around peak.We describe the radiative transfer models of SNe Ia in the nebular phase, which were first presented by Shingles et al. (2022), in Section 3.3.

Fitting the NIR data
To determine whether a light curve displays a plateau, and if so, when the transition onto the plateau occurs, we performed one-and two-component fits to the light curves between 30 -500 d using a Markov Chain Monte Carlo (MCMC) using the package (Foreman-Mackey et al. 2013).If a two-component fit is preferred over a one-component fit, we classify the light curve as having a plateau.We opted to use MCMC to perform these fits to obtain robust estimates of the uncertainties on each parameter.
For the one-component fit, we fit the following equation to   (), the magnitude in filter  at time : where  2 is the slope and  2 is the y-intercept.For the two-component fit, we implemented the same method as that used by Anderson et al. (2014) for characterising the light curves of SNe II.The twocomponent fit is described by the following piece-wise function: Where  2 and  2 are the same as for the one-component fit, and  1 and  1 are the slope and y-intercept of the function prior to the transition onto the plateau.The time of the onset of the plateau in filter ,   0 is defined as:  15) SOFI * We found a very large uncertainty on the magnitude (18.5 ± 3.0 mag) for one J-band image of SN 2021aefx at MJD = 59661.0d.This is likely because the source was faint and the exposure time (60 s) was not sufficient.The next data point at MJD = 59816.3d has a similar magnitude (18.89 ± 0.24 mag) but was exposed for 1080 s and has a significantly smaller uncertainty.We quote this data point as an upper limit at 15.5 mag. to ensure that the two linear components intersect at   0 .We ran an MCMC using 10 walkers for 10,000 iterations and uninformative priors.To avoid biasing the estimates of the slope during the plateau, we exclusively used data taken in J/H or F125W/F160W.We required at least four data points to perform the two-component fit since we are fitting for four parameters, and we required at least one data point at < 150 d and one at > 150 d to ensure we are sampling the phase ranges at either side of the expected transition onto the plateau.Only SNe 2001el, 2011fe, 2012ht, 2014J, 2018gv, 2021pit, and 2021aefx had sufficient data coverage to perform both one-and two-component fits across the range 30 -500 d.For the rest of the sample, there is not enough data to determine the plateau onset and we only performed one-component fits between 150 -500 d to find a single slope and y-intercept ( 2 ,  2 ).For the objects where only the one-component fit was possible, we limited the phase range to 150 -500 d because it is unclear whether these SNe display a plateau phase, and we want to ensure we do not include data before the transition onto the plateau.At least two data points were required per band per SN to perform the one-component fit.
On short timescales, the photometric uncertainty dominates over the temporal evolution, which results in highly uncertain estimates of the slope.We therefore required at least two data points to sepa- We present an overview of the J-, H-, and K  -band light curves in absolute magnitude for all the SNe Ia in our sample.We note that the NIRI photometry of SN 2020uxz is taken in K rather than K  .We also include a comparison to the sub-M ch model with 8x heatboost that best matched SN 2013ct from Shingles et al. (2022).The model is scaled to our photometry in the J band.rated by at least 25 d, which reduced the number of SNe Ia with a measurement of the decline rate to 14.The minimum spacing of 25 d was determined by comparing the expected evolution with the expectation fluctuation within uncertainties.The mean uncertainty on the magnitude across our sample is 0.3 mag.The decline rate in the K  band measured across the whole sample is 1.2 ± 0.2 mag / 100 d, meaning that a change of ∼0.3 mag would be expected to occur across approximately 25 d.
We used reduced- 2 ( 2 red ) to describe the quality of a fit and we used the Akaike Information Criterion (AIC; Burnham & Anderson 2004) to determine the rank of the one-and two-component fits.AIC penalizes extra degrees of freedom to avoid over-fitting the data, and is defined as follows: where  is the likelihood and  is the number of free parameters.For comparing the one-and two-component models, we have  = 2, 4, respectively.Therefore, the two-component model is penalised for its two additional degrees of freedom by four AIC units.If the one-and two-component models differ by more than 2 AIC units, the model with a lower score was deemed the better fit.The best matching value is taken from the 50th percentile, and the uncertainties were taken from the 16th and 84th percentiles of the marginalised distribution.We used the best matching fits to determine the properties of the plateau for each SN.The decline rate during the plateau phase is taken as the slope ( 2 ), which we quote in units of mag / 100 d.The average magnitude is calculated from  2 and  2 of the best matching fit.
To minimise the impact from poorly sampled light curves, we also performed the one-and two-component fits for the full combined sample in absolute magnitude for each filter.Since there is minimal intrinsic scatter in the NIR, this should give a good estimate of the average decline rate for each filter.

SALT3 light curve fits
In order to determine the general properties of each SN, we fitted optical light curves taken from the literature with the package , using the SALT3 model (Kenworthy et al. 2021).The sources of the optical light curves are listed in Table A1.For SNe with no published optical data, we used preliminary photometry from LCO provided by the Global Supernova Project (GSP).We excluded any UV or NIR data because SALT3 is not well trained at those wavelengths, and we restricted the data to between −10 d to +40 d.The SALT3 parameters derived from these fits ( 1 , a metric of the light curve stretch, and , a measure of the colour at peak) are presented in Table A1.There was no optical light curve available for SN 2013ct so we were not able to derive SALT3 parameters.

Comparison to radiative transfer models
We compared our sample to the sub-M ch SN Ia models of Shingles et al. (2022).These models use the Shen et al. (2018) model of a detonation of a 1  WD and evolve the post-explosion composition using the radiative transfer code .Earlier models by Fransson et al. (1996) predicted a strong decline in the optical as flux is redistributed to the NIR, which was not matched by observations.The improved treatment of non-local scattering and fluorescence by Fransson & Jerkstrand (2015) alleviated some of the discrepancies between the models and the observations, but no light curves were published for direct comparison.The models of Shingles et al. (2022) use a modified treatment of non-thermal energy deposition in which the the energy loss to free electrons is artificially boosted as a way to lower the ionisation state.With this modification, the models are able to reconcile the strength of the [Fe ] features, which are generally under-produced by sub-M ch models.Others have suggested that clumping of the ejecta is required to reduce the ionization state (Wilk et al. 2018).The sub-M ch model with a plasma loss rate increased eight-fold (model 1) is best able to reproduce the nebular NIR spectrum of the normal SN 2013ct (see fig. 6 in Shingles et al. 2022).In this work, we present a time-extended version of the sub-M ch -heatboost8 model.For further details of the the model, we refer the reader to Shingles et al. (2022).We also include the other three sub-M ch models presented by Shingles et al. (2022) (sub-M ch -heatboost×4, sub-M ch , sub-M ch -AxelrodNT), referred to from hereon as models 2, 3, and 4, respectively.However, because SN 2013ct is a normal SN Ia and is included in our sample, we focus on the best-matching model to its nebular spectrum (model 1).In Fig. 1 the model is scaled to the J-band photometry from our sample.

RESULTS
We constrain the onset of the plateau for a sub-set of SNe Ia in Section 4.1.1.In Sections 4.1.2and 4.1.3,we present the decline rates and average magnitudes during the plateau of the SNe Ia in our sample.We compare the NIR plateau properties to SN properties at peak in Section 4.2.

Properties of the NIR plateau
We analyse the photometry presented in Section 2 and shown in Fig. 1 using three metrics: the onset of the plateau, the decline rate, and the average magnitude during the plateau.These metrics are derived from the fits either to each individual SN or to the sample as a whole in each filter, as described in Section 3.1.2001el, 2011fe, 2012ht, 2014J, 2018gv, 2021pit, and 2021aefx have data before and during the plateau, enabling us to constrain the phase of the onset of the plateau.The onset of the plateau is calculated by fitting a two-component linear fit to the light curves, as described in Section 3.1 and shown in Fig. 2. We also fit each SN with a one-component fit and compare the result to the two-component fit using the AIC.If the one-and two-component fits differ by more than two AIC units, the model with the lower AIC is deemed significantly better.

SNe
SN 2021pit, which is the best sampled SN along the transition onto the plateau, is best fit with two components in J and H, and yields   0 = 130 +30 −70 d and   0 = 160 +30 −40 d.The K  band is best fit with a single-component decline (ΔAIC = 5).The uncertainty on the time of the onset of the plateau is large and we are unable to constrain the onset to the order of a few days, likely due to the gradual nature of the transition.
The results from the fits to SNe 2011fe, 2012ht, 2014J, 2018gv, and 2021aefx are summarised in Table 4.All Jand H-band light curves are best matched by a two-component model.Only SNe 2001el and 2021pit have sufficient K  -band data to perform one-and twocomponent fits, and they were both best matched by a one-component model suggesting no plateau exists in this band.
We find that the   0 and   0 values for the other SNe are consistent with those derived for SN 2021pit.We note that SN 2012ht has pretransition data in F160W whereas the post-transition data is taken with H, therefore the parameters derived describing the transition onto the plateau should be treated with caution.Similarly, SN 2021pit has post-transition data from HST in F160W.However, repeating the fit excluding the HST data produces consistent results (  0 = 140 ± 30 d).
We calculate the weighted mean of  0 for all measurements across one filter, taking into account the uncertainties, and find   0 = 90 ± 20 d and   0 = 130 ± 20 d, where the uncertainties are quoted as the standard deviations.This implies   0 and   0 are consistent.This disagrees with the trend for the secondary maximum, where the second peak occurs in H before it occurs in J (Kasen 2006, Dhawan et al. 2015).However, a larger, better sampled collection of SNe is required to reduce the uncertainties on the time of the transition and test if the time of transition is truly consistent between the J and H bands.
To increase the sample size, we repeat the same one-and twocomponent fits for the full combined sample (Fig. 2).We find that the Jand H-band data are best fit with two components (ΔAIC = 1137 and 3832), with   0 = 120 ± 10 d and   0 = 140 ± 10 d.The K  band is best fit with a single, constantly declining component (ΔAIC = 580).

Decline rate during the plateau
In Fig. 2 we show two-component linear fits, fitted to all the SNe Ia simultaneously in absolute magnitude in each filter.By fitting the full combined sample, the influence from a single, potentially poorly   4):  1 and  2 are the slopes prior to and during the plateau.We also include  2 for SNe Ia that exclusively have data during the plateau, for which we performed only one-component fits.Column (5):  0 is the phase at which the SN transitions onto the plateau.Column ( 6) & (7):  2 red values of the one-and two-component fits, describing the quality of the best-matching fit.Columns (8): The difference between the AIC values for the one-and two-component fits (ΔAIC = AIC one−comp.− AIC two−comp.).If the AIC values of two models differ by more than 2 units, the model with the lower AIC value is deemed significantly better, and the parameters for that model are quoted.* The pre-and post-transition data have contributions from both space and ground based telescopes.sampled SN, is minimised.In the J and H bands, the decline rates of the second component (during the plateau) are   2 = 0.4 ± 0.1 and   2 = 0.5 ± 0.1 mag / 100 d, respectively.These are inconsistent with zero at a >3 confidence level, meaning that the decline does not cease completely during the plateau.However, by comparing to the decline rate prior to the plateau (  1 = 4.3 ± 0.6 and   1 = 3.6 +0.4 −0.3 mag / 100 d) it is clear that the decline slows significantly.The decline rate in K  is inconsistent with zero at a >6 confidence level, since this band is best fit by a single component with a continuous decline (  2 = 1.2 ± 0.2 mag / 100 d).
In Fig. 3 we show the decline rate for each individual SN as a function of the absolute B-band magnitude at peak ( B max ) and the mean phase of the observations.In the J band, most SNe Ia have a slope consistent with zero.The decline rate averaged across the SNe Ia in the J band is 0.4 ± 0.4 mag / 100 d, where the uncertainty is the standard deviation weighted by the individual uncertainties.In the H band, the average decline rate is 0.5 ± 0.3 mag / 100 d.The K  band behaves differently from the other two bands, with an average decline rate of 1.3 mag / 100 d and a weighted standard deviation of 0.3 mag / 100 d.
When comparing the decline rates between J and H for each SN Ia in the sample, they are consistent for six SNe Ia.One SN Ia has a steeper decline in J, whilst two have a steeper decline in H.Those with a steeper decline in the H band have observations limited to < 250 d, whereas those with a consistent decline rate, or shallower decline rate in H, have observations taken later than 250 d.This points to an evolution in the decline rate in the H band across the plateau, with a steeper intial decline in H which levels off with time.This evolution was also apparent in the well-sampled light curves of SN 2017erp and SN 2018gv (Graur et al. 2020), where at the start of the plateau phase the decline rate decreases, but near the end of the plateau phase the decline rate begins to rise again.Unfortunately, no well-sampled J-band light curve is available for comparison, but we refer the reader to Appendix B for an analysis of this evolution for model 1.

Average magnitude during the plateau
In Fig. 4 we show the average magnitude as a function of the light curve stretch ( 1 ), as well as the mean phase during which the data were taken.SNe Ia with only a single data point during the plateau are included but are indicated by markers without a black outline.For the J and H bands a single data point should give a reliable estimate of the average magnitude during the plateau due to the approximately flat decline rate in these two filters.However, the K  band estimates for these SNe Ia are more uncertain due to the steeper decline in this filter.

Correlations between plateau properties and SN properties at peak
Graur et al. ( 2020) find that the average magnitude during the plateau in the H band scales with Δ 100 () (the decrease in magnitude between peak and 100 d after peak in the H band) and Δ 15 ().In the following section we explore the correlations between the plateau properties,  B max , , and  1 (available in Table A1).To measure how strongly two variables are linearly related, we use Pearson's correlation coefficient, .The significance of the correlation is measured by the -value, with <0.05 indicating a statistically significant correlation.We find a significant correlation between  1 and the average Jand H-band magnitudes during the plateau ( = −0.54,−0.55 and  = 0.038 and 0.036, respectively), implying that SNe Ia with broader light curves (larger  1 values) are intrinsically brighter in J and H during the plateau (Fig. 4).This trend agrees with the correlation found by Graur et al. (2020) for the H band.We find no statistically significant correlation between the average magnitude during the plateau and  1 in the K  band.
The average magnitude during the plateau is driven predominantly by the slope of the decline prior to the transition onto the plateau, a metric that can be approximated by Δ 100 (), as shown by Graur et al. (2020).Δ 100 () shows a weak correlation with Δ 15 (), as shown in fig. 3 of Graur et al. (2020).Combining these results from Graur et al. (2020) and this work, we suggest that broader SNe Ia (larger  1 , smaller Δ 15 ()) tend to decline less in the period 100 days after maximum in H and therefore have a higher average magnitude during the plateau phase.
One potential source of bias worth considering is that the likelihood of being able to observe a SN Ia during the plateau is a function of its brightness on the plateau.If a SN Ia has a shallower decline after maximum (smaller Δ 100 ()), it will remain brighter during the plateau.Therefore, it is likely that studies of SNe Ia on the plateau are inherently biased and tend to sample the SNe Ia that are brighter during the plateau and lie at the lower end of the Δ 100 () population.
We find no significant correlations between the slope during the plateau and  B max ,  1 , or .However, we note that SNe 2000cx and 2013aa are clear outliers in  B max vs. slope (Fig. 3).We check for a linear correlation between the slope and  B max excluding these two SNe Ia.The result from this fit is shown as the dashed line in Fig. 3.This correlation is significant in H with Pearson's  coefficient = −0.83and -value = 0.002, implying that SNe Ia that are more luminous at peak tend to decline faster during the plateau phase.It is unclear why SNe 2000cx and 2013aa do not follow this trend, but Both the J and H bands have a non-zero decline rate but are consistent with zero within 2, and are best fit with two components.The decline rate in the K  band is inconsistent with zero at a >6 confidence level, and is best matched by a one-component fit.We also show the best matching fits to the SNe Ia with sufficient data (markers are the same as in Fig. 1).both are very luminous at peak ( B max < -19.5).We discuss these objects in more detail in Section 5.3.
The timing of the secondary maximum of SNe Ia shows a strong correlation with the stretch of the light curve (Dhawan et al. 2015, Papadogiannakis et al. 2019), with narrow, fast evolving SNe Ia having an earlier secondary maximum.An increase in the total mass of 56 Ni (corresponding to a smaller Δ 15 ()) delays the onset of the secondary maximum due to the higher temperature of the ejecta (Kasen 2006).We suggest that the NIR plateau is caused by a similar mechanism as the secondary maximum, and we expect similar correlations to hold for the NIR plateau.The timing of the onset of the plateau could therefore be expected to correlate with the stretch of the light curve.We test whether there is any correlation between  1 and the  0 values calculated in Section 4.1.1 and find no statistically significant correlations (-values = 0.8 and 0.1 for the J and H bands, respectively).However, for most SNe the phase of the onset is very poorly constrained due to poor sampling, and we cannot rule out a possible correlation between these parameters.Future studies of SNe Ia with higher cadence observations (< 20 d) around the transition phase (70 -150 d) will help to answer this question.We include measurements taken from a single data point (shown as markers without a black outline), which should give a reasonable estimate of the magnitude in the flatter J and H bands, but should be interpreted with caution for the K  band due to its steeper evolution throughout the plateau.We perform a linear fit (dashed line) in each filter and include the calculated Pearson's coefficient and the corresponding -value.We find a significant (  < 0.05) trend of the average magnitude during the plateau with  1 in the J and H bands.

DISCUSSION
In Section 5.1, we provide a theoretical discussion of the NIR spectral evolution.We then answer the questions raised by Graur et al. (2020): "Is there a plateau in the K  band?" and "Does the H-band plateau consist of two distinct branches?",by analysing the results presented in Section 4. We discuss how the models presented in Section 3.3 compare to our observations in Section 5.2.In Section 5.3, we discuss the peculiar SN Ia sub-types present in the sample and compare their behaviour on the plateau to the normal SNe Ia.

Relating the photometric evolution to spectral features
The NIR spectrum during the plateau contains many forbidden iron group lines.We show the spectral evolution of SN 2014J throughout the plateau in Fig. 5, with the main spectral features indicated (these spectra were previously published by Dhawan et al. 2018 andDiamond et al. 2018).The strength of the lines at 1.54 µm, 1.74 µm (H band), 2.02 µm, 2.15 µm, 2.22 µm, and 2.35 µm (K  band) decrease with time.These lines contain emission features coming from [Co ], [Co ], [Fe ], and [Fe ], although from Fig. 5 alone it is not possible to say which emission lines from which elements dominate each feature.
To learn more about the individual contributions to each emission feature, we use information about the transition probability of each line from the atomic data made available by the National Institute of Standards and Technology (NIST). 4The transitions are optically  thin, so their fluxes are proportional to their upper level population times their emission probability.The upper level populations will be similar if their excitation energies are similar, they have similar statistical weights (g=2J+1, where g is the statistical weight and J is the quantum number representing the combined total angular momentum of the electron), and they are both metastable states (only forbidden downward transitions).Therefore, we use the ratios of transition probabilities as a proxy for line strength ratios if the emission lines originate from the same species and have similar excitation energies for the upper level (Jerkstrand et al. 2015).
The two [Co ] features in the J band (1.27 µm and 1.31 µm) and the two [Co ] features in the H band (1.54 µm and 1.74 µm) all have similar upper energy levels (23 060.95, 23 060.95 , 23 435.93, and 22 721.42 cm −1 , respectively).The 1.27 µm and 1.31 µm lines originate from the same  4  multiplet, but come from states with J = 5/2 and J = 3/2, respectively, meaning that the 1.27 µm feature is expected to be about (5+1)/(3+1) = 1.5 times stronger.The 1.54 µm and 1.74 µm lines come from the same upper state ( 2  9/2 ), so the  ki ratio provides a reliable estimate of the flux ratio of these two lines.The emission line at 1.54 µm has the highest transition probability ( ki = 1.3 × 10 −1 s −1 ), whereas the lines in the J band have transition probabilities of 5.4 × 10 −3 and 2.5 × 10 −3 s −1 , respectively.The line at 1.74 µm has a transition probability of 4.2 × 10 −2 s −1 .The dominant [Co ] features therefore sit in the H band, and this band will be most impacted by the decay from 56 Co → 56 Fe.
The features in the J band show only limited decay with time relative to the H band, which aligns with the lower transition probabilities of the [Co ] features at these wavelengths.This is further supported by Fig. 6, which deconstructs the spectrum of model 1 to show the contributions from different species.This is likely not the only correct model for all SNe Ia, but it demonstrates the commonly identified features whilst also being able to model the plateau behaviour.The model suggests that the J band is dominated by [Fe ], whereas the H band has significant contribution from [Co ].The feature at 1.74 µm is composed of three emission lines from [Fe ], [Fe ], and [Co ].The feature at 1.74 µm is dominated by [Fe ], but the model demonstrates that it also has significant contribution from [Co ], explaining its decay with time.This is in agreement with previously identified features in SNe Ia.The K  band is dominated by [Fe ] features (see Fig. 6).

5.1.1
Is there a plateau in the K  band?Graur et al. (2020) speculated that based on the synthetic photometry of SN 2014J, the NIR plateau does not extend to the K  band.In Section 4, we presented additional data in the K  band supporting this conclusion.Here, we rationalise the lack of a plateau in the K  band by referencing the spectroscopic evolution of SN 2014J as a representative of a normal SN Ia (see Fig. 5).
As shown in Figs. 5 and 6, the -band is dominated by an [Fe ] We note that the photometry in fig. 2 in Graur et al. (2020) is scaled to   max whereas the photometry in Fig. 1 is not scaled because H-band data around peak is not available for all the SNe Ia in our sample.A direct comparison between the plots is therefore not possible, but we note that the H-band data in this paper is not separated into two different branches.This could mean that the magnitudes of the plateau in the H band make up a continuous distribution, but only the extremes of this population were sampled by Graur et al. (2020).
We test this first explanation by comparing the average H-band magnitudes on the plateau of the SNe Ia presented by Graur et al. (2020) to the additional SNe Ia presented in this paper (Fig. 7).The gap between −11.5 and −12.5 mag found by Graur et al. (2020) is populated by the SNe Ia presented in this paper, suggesting that the Hband magnitudes on the plateau represent a continuous distribution.A simple Kolmogorov-Smirnov (KS) test enables us to test whether the two samples are likely sampled from the same distribution.We find a KS-value = 0.4 and -value = 0.5, suggesting that the two samples are most likely drawn from the same population.
Despite not seeing two separate branches in the H-band plateau, the behaviour in the H band is less homogeneous than in the J band.Fig. 5 shows a strong [Fe ]/[Co ]/[Co ] complex present in the H band, which is dominated by [Fe ] and [Co ].This feature decreases in strength throughout the plateau due to the continued decay of 56 Co to 56 Fe (Childress et al. 2015, Flörs et al. 2018), which would suggest that there should be some decrease in flux in the H band during the plateau phase.This aligns well with the model predictions, which suggest that the J band decays slower during the plateau phase than the H band.We would expect the decline rate in the H band to flatten with time, as the relative contribution of the [Co ] feature decreases and the decay of these features will have a smaller overall impact on the integrated flux across the filter.Generally, the SNe Ia with observations taken at later phases have shallower declines (with the exception of SN 2020uxz), suggesting that the variation seen in the H band could be driven by whether the observations are taken during the early or late stages of the plateau.

Comparing models to the observations
The light curve of model 1 is shown in Fig. 1, where it is scaled to the J-band photometry of our sample, the wavelength range best matched by the model presented by Shingles et al. (2022).In the H and K  bands, the model under-predicts the magnitude.In the H band the discrepancy is greatest near the beginning of the plateau but lessens with time, whereas in the K  band the offset remains constant.This mismatch between the relative model flux and observed flux in each band is likely due to specific spectral features not being reproduced as well by the models.Figure 6, which is an extended version of fig. 5 from Shingles et al. (2022), highlights that the J band is dominated by an [Fe ] complex spanning 1.22 -1.36 µm.The H band contains a complex of [Fe ], [Fe ], [Co ], and [Co ].Fig. 5 in Shingles et al. (2022), which compares the model spectrum to the spectrum of SN 2013ct, demonstrates that the model is able to reproduce most spectral features across the J and H bands.However, the feature at 1.54 µm, which is dominated by [Fe ] in Fig. 6, is underestimated.This feature also has contributions from [Fe ], [Co ], and [Co ], which may be underestimated by model 1.
In Fig. 3 we show the average decline rate of model 1 on the plateau, between 150-500 d.We find ∼0.2, 0.3 and 1.1 mag / 100 d in the J, H, and K  bands, respectively.This is in agreement with the observational data regarding the presence of the plateau in the H and J bands, as well as the lack of a plateau in the K  band.Moreover, the decline rate predicted by model 1 sits in the parameter space defined by our sample.Fig. 3 also shows the average decline rates for the other three sub- ch models.All four models fall within the parameter space set by the observed SNe Ia, although in the J and H bands, models 2, 3, and 4 tend to predict steeper declines than the majority of our sample (with the exception of SN 2021wuf).We provide a more detailed analysis of the magnitude evolution of model 1 in Appendix B, including an analysis of the first and second derivatives, to characterise the evolution of the slope as well as the inflection points.

Peculiar SN Ia sub-types on the plateau
The majority of the SNe Ia presented in this paper are classified as "normal" SNe Ia based on their maximum-light spectra (see Table A1), although there are a few exceptions.SN 2021wuf is classified as a 91T-like SN Ia, a subclass that follows the width-luminosity relation (Rust 1974, Pskovskii 1977, Phillips et al. 1993) and is used for cosmology, but with light curves that are generally brighter and slower evolving than normal SNe Ia.They show a preference for exploding in late-type galaxies (Taubenberger 2017).SN 2000cx is a peculiar SN Ia, with properties similar to the 91T-like sub-class but with an asymmetric B-band light curve and a peculiar colour evolution (Li et al. 2001).SNe 2004eo and 2012ht are classified as transitional objects between normal and sub-luminous SNe Ia (Yamanaka et al. 2014).Whether these transitional SNe Ia should be used for cosmology is an on-going debate (Gall et al. 2018, Burns et al. 2018, Dhawan et al. 2022, Harvey et al. subm.).SNe 2000cx, 2012ht, and 2021wuf all have consistent average magnitudes in H during the plateau and follow the trend that narrower SNe Ia tend to be fainter during the plateau.SN 2004eo sits above this trend in H, being more luminous than expected for its measured  1 .However, in J SN 2004eo is consistent with this trend.
We found a correlation between  B max and the slope, and whilst the low decline rate of SN 2012ht fits into this trend, SN 2000cx is a clear outlier (Fig. 3).SN 2013aa, whilst being classified as a normal SN Ia, is also exceptionally luminous and similarly falls outside this correlation.On the other hand, SN 2021wuf shows a steep decline during the plateau, as expected from the correlation, although we note that this measurement is based on only two data points separated by 26 d (a minimum of 25 d is required to calculate a reliable slope).SN 2004eo only has two data points in J, and these are not sufficiently spaced to calculate a decline rate.
A larger sample is required to investigate these trends, but if overluminous SNe Ia tend to have a flatter plateau in J and H, this may imply that there is an additional spectral contribution at these wavelengths supporting their luminosity for a longer period.
It has been suggested that the single-degenerate scenario with a near-M ch WD could be solely responsible for over-luminous 91T-like SNe Ia rather than the normal SN Ia population (Fisher & Jumper 2015, Byrohl et al. 2019, Childress et al. 2015).Previous studies have also found that over-luminous 91T-like SNe Ia show flux excesses at a higher rate than normal SNe Ia, which could point towards interaction with a non-degenerate companion in the single-degenerate scenario (Jiang et al. 2018, Deckers et al. 2022 but see Burke et al. 2022 for an alternative view).

CONCLUSIONS
We present NIR photometry of 24 SNe Ia during the plateau phase.From this extensive data set we are able to measure the average magnitude and slope of the plateau in J, H, and K  .We compare these plateau properties to the properties at maximum light and find a significant correlation between  1 and the magnitude of the plateau in J and H, as well as between  B max and the slope in H. From these correlations we conclude that the main driving factor for the magnitude of the plateau is the luminosity at maximum light, which in turn correlates with the decline in magnitude in H 100 d after maximum, (Δ 100 ()).SNe Ia which are more luminous at peak appear to decline faster during the plateau, although there are clear outliers to this trend.Specifically, the over-luminous SNe in our sample behave differently from the normal SNe Ia during the plateau.Over-luminous SNe Ia appear to decline slower than predicted by the trend found between  B max and the slope, which could imply that there is an additional spectral contribution during the plateau.
We constrain the onset of the plateau to 70 -150 d.The secondary maximum occurs in H before it occurs in J (Kasen 2006, Dhawan et al. 2015), but due to the large uncertainties in our estimates of the transition phase we are unable to determine if this is the case for the plateau.We expect a correlation to exist between the time of the onset of the plateau and the peak luminosity of a SN Ia, akin to the correlation found for the secondary maximum, but this could not be confirmed for our sample.
We compare our photometry to models produced by Shingles et al. (2022) and find good agreement regarding the evolution during the plateau, albeit the models under-predict the luminosity in H and K  .However, the best-matching model has reduced non-thermal ionisation rates which leads to lower ionisation states, but no physical justification for reducing these rates has yet been proposed.
An analysis of six spectra of SN 2014J taken throughout the plateau enables us to explain the presence of the plateau in J and H, as well as the absence of the plateau in K  .The dominant [Fe ] features which remain constant throughout the plateau sit in the J and H bands, whilst the K  band hosts mainly [Fe ] features, which recombine to [Fe ] during the plateau phase.
A very limited number of SNe Ia have NIR coverage during the onset of the plateau.Extending this parameter space by obtaining higher cadence observations (< 20 d) around the transition phase (70 -150 d) will enable us to test whether the timing of the plateau correlates with the magnitude at peak, as is the case for the secondary maximum, although we note that this is often difficult due to visibility constraints from the ground.We strongly encourage follow up of over-luminous SNe Ia to test whether they all decline faster during the plateau than expected, since this might imply these events have a different origin.Finally, obtaining more UV photometry coeval with NIR photometry would enable us to determine if flux truly is being redistributed from the UV to the NIR.

Figure 1 .
Figure 1.We present an overview of the J-, H-, and K  -band light curves in absolute magnitude for all the SNe Ia in our sample.We note that the NIRI photometry of SN 2020uxz is taken in K rather than K  .We also include a comparison to the sub-M ch model with 8x heatboost that best matched SN 2013ct fromShingles et al. (2022).The model is scaled to our photometry in the J band.

Figure 2 .
Figure 2. The observed light curves and two-component MCMC fits for the J (top), H (middle), and K  (bottom) filters.The best fits (one-or twocomponent) to the NIR light curves of the whole sample are shown as black solid lines, with the various MCMC iterations shown as faded grey lines.Both the J and H bands have a non-zero decline rate but are consistent with zero within 2, and are best fit with two components.The decline rate in the K  band is inconsistent with zero at a >6 confidence level, and is best matched by a one-component fit.We also show the best matching fits to the SNe Ia with sufficient data (markers are the same as in Fig.1).

Figure 3 .
Figure 3.The measured decline rate as a function of  B max for the J (top), H (middle), and K  (bottom) filters.The marker colour represents the mean phase of the observations.We include the predictions from the sub-M ch models presented by Shingles et al. (2022), using -19.2 as  B max from Shen et al. (2018).The dashed line shows a linear fit to the sample, excluding SNe 2013aa and 2000cx, which are considered as outliers and are highlighted by a thicker marker edge.We also show the Pearson  -coefficient and the corresponding -values.A significant correlation (  < 0.05) is identified in the H band.

Figure 4 .
Figure 4.The average magnitude measured for each SN Ia during the plateau as a function of  1 , with the colour indicating the mean phase of the observations, for the J (top), H (middle), and K  (bottom) bands.We include measurements taken from a single data point (shown as markers without a black outline), which should give a reasonable estimate of the magnitude in the flatter J and H bands, but should be interpreted with caution for the K  band due to its steeper evolution throughout the plateau.We perform a linear fit (dashed line) in each filter and include the calculated Pearson's coefficient and the corresponding -value.We find a significant (  < 0.05) trend of the average magnitude during the plateau with  1 in the J and H bands.

Figure 5 .
Figure 5. NIR spectral series of SN 2014J, covering 307 -478 d.Left: The wavelength range covering the J and H bands. Right: The K  -band spectrum (note the different range on the y-axis).The key spectral features are marked, and the telluric regions are indicated by grey regions.The flux has been normalised to the [Fe ] feature at 1.65 µm.The grey solid curves show the 2MASS J, H, and K  filter response functions.

Figure 6 .
Figure 6. Figure adapted from Shingles et al. (2022) with showing the contribution from each species to the sub-M ch heatboost8 spectrum at 247 d.Most of the J and H bands are dominated by [Fe ], with some [Co ] contributions in the H band.The K  band is dominated by [Fe ] emission.

Figure 7 .
Figure 7.A histogram comparing the average magnitude of the plateau in the H band between Graur et al. (2020) and this paper.Graur et al. (2020) noticed bimodal behaviour in the average magnitudes of the H band, but the additional data presented in this paper suggests the H-band magnitudes represent a continuous distribution.

Figure B1 .
Figure B1.Top panels: The light curves in J, H, and K  bands of model 1, as well as SNe 2017erp and 2018gv.All light curves are fit using a univariate spline, and the resulting fits are plotted as a dashed line for the model and solid lines for the SNe Ia.Middle panels: The first derivative of the apparent magnitude with respect to time.Bottom panels: The second derivative of the apparent magnitude with respect to time.We denote where the second derivative equals zero with a grey line, since where this line meets the second derivative indicates the inflection point in the light curve.Model 1 and SN 2017erp both reach an inflection point at 317 d, whilst SN 2018gv reaches an inflection point at 327 d.

Table 1 .
Overview of NIR photometry obtained with HST.A machine readable version of this table is available in the online material.Phase is defined as the time since maximum light in the B band (MJD - 0 ). a

Table 2 .
Overview of NIR photometry obtained with Gemini, the New Technology Telescope, and Calar Alto observatory.A full, machine readable version of this table is available in the online material.

Table 3 .
Overview of synthetic NIR Photometry obtained with XShooter.

Table 4 .
Summary of parameters for the fits performed for SNe Ia with data available before and after the transition onto the plateau.We perform both a one-and two-component fit to each light curve.SNe Ia with a maximum separation of less than 25 d between data points are not included (see Section 4.1.2).

Table A1 .
A table summarising the SNe Ia present in our sample, showing the classification of the SN, redshift, distance modulus, and time of maximum light.We also include