Sub-Chandrasekhar progenitors favoured for Type Ia supernovae: evidence from late-time spectroscopy

A non-local-thermodynamic-equilibrium level population model of the ﬁrst and second ionization stages of iron, nickel, and cobalt is used to ﬁt a sample of XShooter optical + near-infrared (NIR) spectra of Type Ia supernovae (SNe Ia). From the ratio of the NIR lines to the optical lines limits can be placed on the temperature and density of the emission region. We ﬁnd a similar evolution of these parameters across our sample. Using the evolution of the Fe II 12 570 –7155 Å line as a prior in ﬁts of spectra covering only the optical wavelengths we show that the 7200 Å feature is fully explained by [Fe II ] and [Ni II ] alone. This approach allows us to determine the abundance of Ni II /Fe II for a large sample of 130 optical spectra of 58 SNe Ia with uncertainties small enough to distinguish between Chandrasekhar mass ( M Ch ) and sub-Chandrasekhar mass (sub- M Ch ) explosion models. We conclude that the majority (85 per cent) of normal SNe Ia have a Ni/Fe abundance that is in agreement with predictions of sub- M Ch explosion simulations of ∼ Z (cid:3) progenitors. Only a small fraction (11 per cent) of objects in the sample have a Ni/Fe abundance in agreement with M Ch explosion models.


I N T RO D U C T I O N
Type Ia supernovae (SNe Ia) are a remarkably homogeneous class of transients which are thought to originate from the explosion of a white dwarf (WD) star in a binary system. Radioactive 56 Ni produced in the thermonuclear explosion of the electron-degenerate matter (Hoyle & Fowler 1960) powers the light curve (Pankey 1962;Colgate & McKee 1969;Kuchner et al. 1994) for several years. Even though SNe Ia have been used as distance indicators for several decades and significantly contributed to our current understanding of cosmology ( CDM and the accelerated expansion of the universe, Riess et al. 1998;Perlmutter et al. 1999), the precise mechanism that leads to the thermonuclear runaway reactions, as well as the progenitor system, remains elusive. E-mail: afloers@mpa-garching.mpg.de Two channels that can lead to the explosion of a WD as a SN Ia have been extensively discussed in the literature. In Chandrasekharmass (M Ch ) explosions the burning front propagates either as a deflagration (e.g. Gamezo et al. 2003;Fink et al. 2014) or a delayed detonation (e.g. Blinnikov & Khokhlov 1986, 1987Khokhlov 1991;Gamezo, Khokhlov & Oran 2005;Seitenzahl et al. 2013). The explosion is naturally triggered by an increase of the central density as the WD accretes material from its companion and comes close to the Chandrasekhar mass limit (M M Ch ). In the sub-M Ch channel, the central temperature of the primary WD never reaches conditions that are sufficient to ignite carbon. However, an explosion significantly below the M Ch may be triggered through dynamical processes such as mergers (e.g. Pakmor et al. 2010Pakmor et al. , 2013Ruiter et al. 2013), double detonations (e.g. Fink et al. 2010;Woosley & Kasen 2011;Moll & Woosley 2013;Shen et al. 2018) or head-on collisions (e.g. Kushnir et al. 2013). For such systems the burning front propagates as a pure detonation Blondin, Dessart & Hillier 2018).
The search for solutions to the SN Ia progenitor problem has been the focus of many studies. For historical supernova remnants one can search for a surviving companion star which was ejected at velocities of a few hundred km s −1 , though no promising candidates have been found so far (Kerzendorf et al. 2018a for SN 1006, Kerzendorf et al. 2018bfor SN 1572. Non-degenerate donor stars of SNe Ia in nearby galaxies should also be visible in deep images as their brightness increases by a factor of ∼10-10 3 , though again, no donor stars have been found for a sample of the closest SNe Ia in recent times Bloom et al. 2012;Shappee, Kochanek & Stanek 2013).
The growth of a WD star to the M Ch limit requires a steady transfer of material from the companion (Nomoto, Thielemann & Yokoi 1984). Material which was expelled from the companion but not accreted on the WD enriches the circumstellar material (CSM). In few cases evidence for such a CSM has been detected (Hamuy et al. 2003, Deng et al. 2004Harris et al. 2018, Graham et al. 2019Vallely et al. 2019, Kollmeier et al. 2019. However, the bulk of Ias do not exhibit any evidence for CSM interaction. When the blast wave from the SN explosion runs through the CSM, electrons are accelerated to relativistic speeds and produce radio emission through synchrotron radiation (Chevalier 1982(Chevalier , 1998Chevalier & Fransson 2006). For a nearby SN Ia such as SN 2011fe radio emission should be observable if it exploded in the single degenerate channel. However, no radio emission was found by Horesh et al. (2012) for SN 2011fe [but see also Nugent et al. (2011) for a counter argument].
One can also distinguish the two channels from direct observations of the aftermath of the explosion itself. Unfortunately, the uniformity of explosion model predictions of SNe Ia makes this a challenging task. A promising difference between M Ch and sub-M Ch models is the mass fraction of neutronized species produced in the explosion. While the progenitor metallicity affects how much neutron-rich material can be produced in both channels, additional neutrons are only available for explosions close to the M Ch due to the high central densities (ρ cen ∼ 2 × 10 9 g cm −3 ) which allow electron capture reactions to take place (Iwamoto et al. 1999;Seitenzahl et al. 2013).
X-ray spectroscopy of SN remnants in the Milky Way (MW) and the Large and Small Magellanic Clouds (LMC and SMC) allowed Park et al. (2013), Yamaguchi et al. (2014Yamaguchi et al. ( , 2015, Martínez-Rodríguez et al. (2017), and Seitenzahl et al. (2019) to estimate the fraction of the neutron-rich stable iron-peak isotopes 55 Mn and 58 Ni. They find considerable differences across their sample, but the number of objects for which such a study can be done is limited.
In this work, we are interested in the composition of the ironrich ejecta of SNe Ia. Theoretical explosion models contain the following isotopes in the central region: (a) 56 Ni, which is the most abundant radioactive isotope and responsible for the heating of the ejecta. It decays within a few days (t 1/2 = 6.075 d) to 56 Co, which in turn decays (t 1/2 = 77.2 d) to stable 56 Fe. 56 Ni can be produced in NSE (Nuclear Statistical Equilibrium) without an overabundance of neutrons (Y e = 0.5) or high densities (Hoyle & Fowler 1960). In our analysis, we treat 56 Ni as a reference point and give other abundances in fractions of the 56 Ni mass.
(c) Stable 54,56 Fe which is directly synthesized in the explosion and not a daughter product of radioactive decay. M Ch explosions produce M54,56 Fe /M56 Ni > 10 per cent, while most sub-M Ch models have M54,56 Fe /M56 Ni < 10 per cent (see references in 'b').
(d) Stable 58 Ni which is synthesized in the explosion. Sub-M Ch explosions contain M58 Ni /M56 Ni < 6 per cent and M Ch explosions have M58 Ni /M56 Ni between 8 and 12 per cent (see references in 'b').
Contributions of slowly decaying neutron-rich material (e.g. 57 Co -the daughter product of 57 Ni) to the quasi-bolometric light curve of the nearby SN 2011fe at >1000 d after the explosion were investigated by Shappee (2017), Dimitriadis et al. (2017), and Kerzendorf et al. (2017). This method was used for other nearby transients SN 2012cg (Graur et al. 2016), SN 2013aa (Jacobson-Galán et al. 2018, SN 2014J (Yang et al. 2018), SN 2014lp (Graur et al. 2018b), and SN 2015F (Graur et al. 2018a). However, the physical processes relevant at such late phases (e.g. ionization/recombination) have long time constants and their onset is poorly constrained by the data (Fransson & Jerkstrand 2015). In particular, it remains unclear what fraction of the radioactive decay energy is converted into optical photons as the majority of the energy is expected to come out in the mid-IR.
SNe Ia complete their transition into the nebular phase roughly half a year after the explosion when the ejecta become fully transparent to optical and near-infrared (NIR) photons and the bare iron core which gives insight into the explosion physics is visible. Nebular phase spectral models build on the early work of Axelrod (1980) and many authors over the years (Kozma & Fransson 1992;Spyromilio et al. 1992;Kuchner et al. 1994;Kozma et al. 2005;Mazzali et al. 2007;Fransson & Jerkstrand 2015;Botyánszki & Kasen 2017;Diamond et al. 2018;Maguire et al. 2018) with spectral synthesis codes of varying complexity.
The method presented herein enables the use of optically thin Ni and Fe lines at optical and NIR wavelengths to constrain the fraction of neutron rich material in the ejecta. This approach increases the number of objects for which the analysis can be performed by about one order of magnitude compared to the number of optical + NIR spectra currently available. The analysis is made possible by using a small sample of optical + NIR spectra to determine the relative strength of the NIR 12 570 Å to the optical emission 7155 Å lines of Fe II. With this relation we can model optical nebular spectra which do not have NIR observations. From the fits to observed late spectra we determine the stable Ni to Fe ratio. We estimate the systematic uncertainty of the method and show that emission lines from singly ionized iron and nickel are sufficient to model the 7 200 Å emission feature. Finally, we discuss the implications of the determined Ni to Fe ratio of the large sample of more than 100 spectra on the various explosion model predictions.

O B S E RVAT I O N S
We extend the XShooter sample of nearby SNe Ia in the nebular phase from Maguire et al. (2018) with SN 2015F (PI M. Sullivan, program ID 095.A-0316; PI R. Cartier, program IDs 096.D-0829, 097.D-0967, and 098.D-0692) and SN 2017bzc (PI L. Dessart,. The epochs of the additional spectra range from ∼200 to ∼420 d after the explosion. An overview of the observations for these two supernovae is given in Table 1. XShooter is an Echelle spectrograph with three arms (UVB, VIS, and NIR) covering the wavelength range of ∼3 000-25 000 Å located at the Very Large Telescope (VLT). The resolution of the individual arms depends on the slit widths. For the observations presented in this work slit widths of 1.0 (UVB), 0.9 (VIS), and 0.9 (NIR) arcsec have been used. The corresponding resolution of the three arms is therefore 5400, 8900, and 5600, respectively. The spectra were reduced using the ESO pipeline with the XShooter module, producing flux-calibrated one-dimensional spectra in each of the three arms (Modigliani et al. 2010;Freudling et al. 2013). We also used a custom post-processing pipeline to combine the rectified 2D-images, perform the sky-subtraction and extract the spectrum (https://github.com/jselsing/xsh-postproc).
Nebular phase spectra of SNe Ia exhibit a number of broad (≈ 7000 to 9000 km s −1 ) emission features (see Fig. 1 (Spyromilio et al. 2004;Flörs et al. 2018). The emission feature at around 13 000 Å is identified as the 12 570 Å a 6 D-a 4 D multiplet of [Fe II]. The double peaked feature around 16 000 Å is composed of a blend of [Fe II] and [Co II] lines of the multiplets a 4 F-a 4 D and a 5 F-b 3 F, respectively. Redwards of the strong telluric absorption feature at ∼ 18 500 Å we also detect the a 2 F 7/2 -a 4 F 9/2 line of [Ni II] in spectra with high SNR .
In the optical we see blends of singly and doubly ionized Fe, Co, and Ni. The strong feature at 4700 Å originates mainly from the 3d 6 5 D-3 F multiplet of [Fe III] (Axelrod 1980;Kuchner et al. 1994). The broad emission centred around 5900 Å is primarily due to [Co III] in the a 4 F-a 2 G multiplet. The identification of the 5 900 Å [Co III] feature is secured by the fact that the relative strength of this feature with respect to e.g. the [Fe III] 4 700 Å feature decreases with time as predicted by radioactive decay of 56 Co (Kuchner et al. 1994;Dessart et al. 2014;Childress et al. 2015). Near the 7 200 Å region the spectra exhibit emission lines of the [Fe II] multiplets a 4 F-a 2 G and a 6 D-a 4 P and the [Ni II] multiplet z 2 D-a 2 F. The identification of the various emission lines in the optical and NIR of SNe Ia in the nebular phase has been extensively discussed in the literature and is considered secure. A detailed overview of the strongest emission lines is given in Table 2.

Summary
We have determined that the line ratio of the 12 570 to 7155 Å [Fe II] lines in the nebular spectra of SNe Ia evolves with supernova age in a predictable log-linear manner. We assume that evolution is valid for supernovae for which we only have optical coverage. The range of electron densities and temperatures that give rise to a given ratio is determined by the atomic data for these transitions. For each epoch we thus have prior knowledge of the range of n e and T. This range is used to determine the ratio of the emissivity per atom for the 7155 [Fe II] to 7378 Å [Ni II] lines and thus determine the range of mass ratios of nickel to iron based on optical data alone at any given epoch.

The model
We use a one-zone model as described in Flörs et al. (2018). We extend the model to include all first and second ionization stages of iron, nickel and cobalt (see Table 3). For this set of ions we solve the non-local-thermodynamic-equilibrium (NLTE) rate equations and compute level populations. Throughout this work, we redshift and extinguish the spectral models instead of correcting the observed spectra. In Table B1, we show the redshift and reddening applied to our models. For the reddening correction in our models we adopt Cardelli, Clayton & Mathis (1989). The strength of the reddening is strongly constrained by the presence of a number of lines arising from the same upper level in different ions (e.g. 12 570 and 16 440 Å of [Fe II]).
We assume that thermal emission is the dominant source of light from the start of the nebular phase until ∼ 500 d after the explosion. During this phase, the ejecta are transparent for optical and NIR photons, allowing us to ignore radiative transfer effects. We also do not consider non-thermal excitations as the energy going into this channel at the relatively high electron densities we determine is also very low (Fransson & Chevalier 1989). We do not include charge exchange and time-dependent terms in the NLTE rate equations. For the set of ions given in Table 3 we solve the NLTE rate equations to obtain the level populations of the ions, which are used to determine the line emissivities.
We compare our parametrized model M to the XShooter observations D described in Section 2 using the approach from Czekala et al. (2015). The likelihood function contains a correlation matrix C which has the uncertainties of the pixels as diagonal elements and the correlations of nearby pixels on the off-diagonals: To account for systematic imperfections of the model (e.g. line profiles are not Gaussian far from the line centre), we use Gaussian processes with a Matérn kernel to add an additional noise term in the correlation matrix at the location of the feature edges (see Czekala et al. 2015). This prevents the sampling algorithm from only choosing a narrow set of parameter values, which yield a better fit in regions where the model is systematically unable to fit the observations. We employ flat priors for all parameters of the model. The upper and lower bounds of the flat priors are chosen in such a way that the posterior parameter distributions are not truncated. We use nested sampling to find the posterior distributions of the parameters of the model that yield good fits with the observed spectrum (https://github.com/kbarbary/nestle, see also Shaw, Bridges & Hobson 2007). Fig. 1 presents the fits results for the spectra given in Table 1. The red line indicates the mean flux of all fit models at each wavelength while the orange shaded area marks the 68 per cent uncertainty of the fit. Fit results for the previously published spectra of the XShooter sample are shown in Flörs et al. (2018). An exemplary zoom into the fit of SN 2017bzc at +215 d is shown in Fig. 2.
For each spectrum we can use the posterior distribution of the model parameters to compute line emissivities of all lines of singly and doubly ionized Fe, Ni, and Co. In this work, we use line ratios of [Ni II] Flörs et al. (2018). While the NIR nebular spectra are easier to model than the optical spectra, the mass ratio of Co II to Fe II changes with time and the number of spectra with NIR coverage is quite limited. In this work we want to make use of several decades of optical nebular phase spectroscopy to determine the distribution of the Ni/Fe abundance and compare our findings with predictions from explosion models.

Calibration of optical spectra of SNe Ia
To determine the Ni II / Fe II mass ratio we compute the ratio of the 7378 Å [Ni II] and the 7155 Å [Fe II] lines (see Fig. 2 panel c).
The conversion of line emissivities to emitting masses requires knowledge of the temperature and density of the emitting material. The one-zone-model employed in this study does not allow us to disentangle these two parameters. However, we find that the evolution of the ratio of the strongest Fe II line in the NIR (12 570 Å) and optical (7155 Å) is very similar across our sample of optical + NIR spectra (see Fig. 2 panels c and e for these lines). This seems to be a natural evolution from high temperatures and high densities towards lower values. Due to the decreasing temperature it becomes more difficult at late epochs to excite the levels giving rise to optical transitions, thus increasing the ratio of the NIR to optical lines. We fit a simple linear relation through our inferred data points (see Fig. 3). The uncertainties of the individual data points are uncorrelated, thus justifying the use of a simple Chi-Square likelihood where In this equation y and y indicate the inferred values and uncertainties of the Fe II 12 570 to 7155 Å ratio for our sample, m is the slope of the fit curve, b is its intersect, and σ is the intrinsic scatter of the population. We add an intrinsic scatter term to the likelihood function that takes into consideration that our sample consists of many different objects. The uncertainty of the fit is then a combination of the uncertainty of slope and intersect and the intrinsic scatter term. We find for the ratio of Fe II 12 570 to 7155 Å log F 12 570 with an intrinsic scatter of 0.06 dex around the best-fitting curve. The choice of the atomic data has only very weak consequences on the inferred NIR/VIS ratio. Translating the NIR/VIS ratio to temperatures/densities does rely on the atomic data, however. The atomic data used throughout this work is given in Table 3. Alone, the optical spectra of SNe Ia do not allow us to constrain the temperature and density of the emitting material in any meaningful way -we can obtain good fits for a wide range of temperatures and densities. However, we notice that for a given Fe II 12 570 to 7155 Å ratio only specific tracks in the temperature/density space are possible. The inference uncertainty of the NIR/VIS ratio translates into a curve with non-zero width in the temperature/density space. The measurement of the NIR/VIS line ratio is considered robust -no other strong lines are present in the 12 500 Å feature, and in the 7000 Å region only Co III of the iron group elements has a weak contribution. We exclude the extremes in the temperature/density space (see grey-shaded areas in Fig. 4) by fitting the many lines of singly and doubly ionized material at optical and NIR wavelengths. Each of the curves in Fig. 4 corresponds to one value of the NIR/VIS ratio. We can thus determine a range of temperatures and densities of SNe Ia in the nebular phase assuming that the Fe II 12 570 to 7155 Å ratio evolves as the red curve in Fig. 3 with a 1σ uncertainty of 0.06 dex. We add this constraint as a Gaussian prior into the likelihood function of our Bayesian fit model.

Determination of the Ni to Fe ratio
The temperature range is significant and thus the uncertainty in the absolute masses is large. A temperature difference of only a few hundred Kelvin can lead to an emitting mass that is different by a factor of a few. However, a more robust quantity is the mass ratio of ions of the same ionization stage. Under the assumption that the same emitting region gives rise to these lines Maguire et al. 2018) the physical conditions of the ions (temperature and electron density) are similar. Using a mass ratio also negates the effect of the rather unknown distance to the SN host galaxy and significantly reduces the effect of the emitting temperature.  For a given temperature and density we can directly infer the ratio of the number of emitting Fe II and Ni II ions required to match the observed flux ratio of the 7155 and 7378 Å lines (see Fig. 2 panel c) for the lines and Table 4 for the resulting Ni/Fe mass ratios. Temperatures and densities that yield a good fit can be found if a NIR spectrum is available. For spectra that lack this additional information we have to use the relation obtained in Section 3.3. We discuss the additional uncertainties from using the fit relation instead of the full optical + NIR spectrum in Section 4.2.

The Fe II NIR/VIS ratio
In Section 3.3, we derived a relation between two of the strongest Fe II lines that are observed in nebular spectra of SNe Ia. Our extended XShooter sample now contains 14 spectra of 9 different SNe. The ratio of the NIR 12 570 and the 7155 Å lines evolves similarly for all objects in our sample. In physical terms, the ratio  of these lines is a direct measure of the cooling and expanding Ferich ejecta. The relation does not depend on the collision strengths but only on the transition rates of Fe II. These are well known, as can be seen from the match of the fit models and the observed spectra in regions where only Fe II emission is present. Additionally, the Fe II NIR/VIS relation as presented in this work is not just the result of a possible oversimplification of our one-zone model. It is obtained by effectively de-blending the lines of singly and doubly ionized iron,    ratio for a given epoch t i has the same width as the fit curve in Fig. 3. In general, fitting the optical spectrum with the use of the NIR/VIS relation does not necessarily prefer the same ratio as fitting the full optical and NIR spectrum. As a result, we obtain different posteriors for the density and temperature for the two fitting methods. It seems that the optical is more sensitive to different regimes of the electron density and temperature than the combined optical and NIR spectrum. On average, the use of the Fe II 12 570 to 7155 Å fit  relation instead of the NIR spectrum leads to a systematic difference of σ sys = −0.0033 +0.0037 −0.0041 . The use of the NIR/VIS relation therefore results in mostly smaller M Ni / M Fe ratios by about 0.0033 within the 68 per cent confidence interval. We consider this a systematic uncertainty that adds to the statistical uncertainty linearly.

Time evolution of the Ni/Fe ratio
Even though the amount of 58 Ni produced in the explosion is fixed for a single object, the ratio of Ni/Fe changes with time (Fe being the daughter product of 56 Co decay, which at early times has not completely decayed yet). Only after ∼300 d (4 × t 1/2, 56 Co→ 56 Fe ) the Ni/Fe ratio remains almost constant.
For supernovae that have several observations during the nebular phase we can test whether our modelling yields consistent Ni to Fe ratios (i.e. that the slope of the data points follows a single theoretical explosion model prediction). In Fig. 6, we normalize the Ni to Fe to the value at t = t ∞ to make it easier for the reader to see the slope of the measured data points. A flat series of data points indicates that the evolution with time behaves according to the expected yields from the radioactive decay of 56 Ni. For objects with both optical and NIR data the full spectrum is fit while for objects with optical data only the method described herein is used to provide the range and evolution of n e and T.
The evolution of the Ni to Fe mass ratio for objects with multiple observations during the nebular phase is consistent with pure radioactive decay within the statistical uncertainties. A much shallower or steeper slope of the NIR/VIS ratio would lead to nonflat evolutionary curves of the Ni/Fe ratio. The only object that shows an evolution of the scaled Ni/Fe mass ratio is SN 1998bu. Roughly 270 d past its B-band maximum the inferred Ni/Fe mass ratio increases by about 15 per cent and settles on this new value for the remaining observations. Such a behaviour could be the result of a light-echo contribution to the nebular spectrum, as was found for SN 1998bu by Cappellaro et al. (2001).

A possible contribution of Calcium at 7200 Å?
The method presented in Section 3 relies on the assumption that only [Fe II] and [Ni II] contribute to the 7200 Å feature. If emission from another ion (e.g. Ca II]) contributes substantially to this feature, our measurement will be systematically wrong as the true contribution of Ni to the feature is lower than estimated from our model. Some NLTE radiative transfer calculations of SNe Ia in the nebular phase predict a non-negligible flux of Ca II] emission at λλ 7291.5, 7323.9 (Botyánszki & Kasen 2017;Wilk, Hillier & Dessart 2019). If the emitting region is a spherical shell at high velocities outside the iron core, the profile would be flat-topped. Such a plateau of Ca II] emission would raise the overall flux level in the 7200 Å region without changing the characteristic double peaked shape of the feature.
We can test whether there is a contribution from other ions by fixing the strength of the Ni II 7378 Å through the 19 390 Å line. The relative strength of the two lines only depends on the extinction and the ratio of the transition rates, as they originate from the same upper level: The observed strength will depend on the extinction. The observations of SN 2015F, one of the closest SNe in the last decade, can be used to further verify this method. We obtained 5 nebular phase XShooter spectra between +181 and +406 d after B-band maximum. The first four epochs (+181, +225, +239, +266 d after maximum) are of exceptional quality and clearly show the 19 390 Å line. The observation at +406 d has a SNR that is insufficient to detect such a weak line, especially as it lies close to a telluric feature. SN 2017bzc was farther away than SN 2015F, but with an integration time of 10 080 s the [Ni II] 19 390 Å line can be seen in the +215 d spectrum.
An overview of the model fits for each of these spectra is shown in Fig. 7. The 19 000 Å feature has not been used to compute the fits. A significant contribution of Ca II in the optical would lead to a much weaker 19 390 Å line, which is in contradiction to our observations. A weak Ca II contribution cannot be ruled out but its effect on the Ni/Fe mass ratio would be very limited. None of the objects with sufficiently high SNR in the 19 000 Å region require any Ca II]. While it is not impossible that some SNe Ia -transitional objects such as the 86G-like or the faint 91bg-like -exhibit Calcium emission in the 7200 Å region, the feature can be explained by only [Fe II] and [Ni II] for the normal and luminous population of SNe Ia (see also Graham et al. 2017).

M Co/Fe from the extended XShooter sample
The additional observations can be used to further the work described in Flörs et al. (2018). As has been noted by several authors Maguire et al. 2018), the singly ionized lines of Fe, Ni, and Co exhibit the same line shift and width. The same holds true for the two additional SNe with nebular phase XShooter observations presented in this work. It is therefore reasonable to assume that the singly ionized species are co-located in the ejecta and share the physical excitation conditions -temperature and density. Our updated model allows us to directly compute the Co to Fe mass ratio without having to use LTE approximations. The effect, however, is quite limited for the NIR lines in question (< 5 per cent). Fig. 8 displays a comparison of the new observations with the ones from Maguire et al. (2018). We find that the three new objects (SN 2012ht, which was not included in the sample of Flörs et al. 2018) have a M Co / M Fe ratio that is consistent with sub-M Ch explosions. Only the spectrum of SN 2012ht allows us to probe the 57 Ni content in the ejecta as all other spectra are significantly younger than 300 d. For them, the ratio instead is a measure of the fraction of stable iron ( 54,56 Fe) to radioactive iron ( 56 Ni decay products).

M Ni/Fe from archival optical spectra
The evolution of the NIR/VIS lines of Fe II allows us to model nebular spectra that cover only the optical wavelength range. We collected 130 spectra of 58 SNe Ia at epochs > 170 d after B-band maximum that have adequate SNR. A full list of all observations used for this study is given in Table B1. The spectra are modelled as described in Section 3. For SNe which have multiple observations in the nebular phase we combine the inferred mass ratios. We report the inferred scaled Ni/Fe mass ratio in Table B1. No corrections (e.g. fitting optical + NIR spectra versus only optical spectra; Section 4.2.1) have been applied to the inferred values. An overview of all objects (XShooter + archival) in our sample is given in Fig. 9. We find that the majority of SNe exhibit Ni/Fe mass ratios below 0.05.
A similar study was conducted by Maguire et al. (2018) for eight objects in their XShooter sample. The same objects are also included in this work, however, a different method for the determination of the abundance ratio is used. Instead of modelling the full spectrum, Figure 6. Inferred mass ratio of Ni and Fe for supernovae with multiple observations between ∼200 and 500 d after the explosion. Explosion model predictions and inferred data points were scaled to the Ni/Fe abundance at t → ∞ to remove the time dependence of M Fe (t) in order to better illustrate the consistency of the method. We assume a rise time of ∼18 d (Ganeshalingam et al. 2011) to compute the time after explosion. Same colours indicate multiple observations of a supernova. Error bars represent the combined statistical fitting uncertainty and the systematic uncertainty from using the NIR/VIS relation if the spectrum only covers the optical wavelength region up to 10 000 Å. Maguire et al. (2018) restrict themselves to the 7200 Å [Fe II] and [Ni II] dominated region. To convert the ratio of the LTE line fluxes to an abundance ratio of nickel and iron, they use average departure coefficients of a W7 model (Nomoto et al. 1984;Nomoto & Leung 2018) at 330 d from Fransson & Jerkstrand (2015). As this model does not allow a determination of the temperature of the emitting material, Maguire et al. (2018) assume temperatures similar to those of Fransson & Jerkstrand (2015) between 3000 and 8000 K.
The inferred abundance ratio of Ni and Fe from Maguire et al. (2018) and this work deviate by about 1.5σ for the same objects. The differences are mainly due to the placement of the (pseudo-)continuum across the 7200 Å region, leading to a different line ratio of Fe II 7155 Å and Ni II 7378 Å. In this work, we opted for a conservative continuum placement as most of it can be explained by a blend of weak lines of other singly and doubly ionized iron group ions (e.g. [Co III], [Fe III]). The departure coefficients corresponding to the allowed range of temperatures and densities (see Fig. 4) of the emitting material are in good agreement with the ones used by Maguire et al. (2018). The use of the Fe II NIR/VIS relation allows us to better constrain the allowed range of the physical parameters of the singly ionized ejecta, leading to reduced uncertainties compared to Maguire et al. (2018). We want to emphasize that both works make use of the same atomic data for the ions in question.

Implications on the explosion mechanism
The various theoretical explosion models of SNe Ia predict different amounts of neutron rich material. In M Ch explosions the amount of synthesized neutron-rich material is determined by two processes: Carbon simmering and neutron-rich burning: Carbon simmering occurs when a WD accretes slowly towards the M Ch . Densities and temperatures in the centre become high enough to ignite carbon, but no thermonuclear runaway happens due to a large convective core that allows for cooling through escaping neutrinos (Woosley, Wunsch & Kuhlen 2004;Piro & Chang 2008). The burning of carbon leads to mostly 13 N and 23 Na, which can subsequently capture electrons which further increases the neutron excess (Chamulak et al. 2008;Martínez-Rodríguez et al. 2016).
Neutron-rich burning to NSE can shift the equilibrium away from 56 Ni to more neutron-rich isotopes ( 54,56 Fe, 57,58 Ni, 55 Mn) (Iwamoto et al. 1999;Brachwitz et al. 2000). Just before the explosion, the high central density of the progenitor WD leads to neutronization through electron capture in the densest region. Neutron-rich NSE burning is only possible if there is a neutron excess in the NSE burning central region.
In sub-M Ch models such processes are not possible as their progenitors cannot reach the required central density. However, an overabundance of neutrons in a high-metallicity progenitor can still lead to the production of neutron-rich IGE (Timmes, Brown & Truran 2003). The fraction of neutron rich to normal material can cover a wide range of values -from close to zero for Z = Z to that of M Ch explosions at several times solar metallicity (Shen et al. 2018). It remains to be seen whether such extremely high-metallicity progenitors really exist.
We focus on the neutron-rich, stable 58 Ni. The presence of a signature line close to 7378 Å reveals that at least some amount of 58 Ni can be found in all normal SNe Ia observed so far. As shown in Fig. 7 the 7200 Å feature can be explained by a blend of mainly [Fe II] and [Ni II]. In principle there will also be varying amounts of stable iron produced during the explosion, but this contribution to the total iron mass is hard to disentangle from the overwhelming fraction of daughter products of radioactive 56 Co.
In contrast to the artificial W7 model (Nomoto et al. 1984), stateof-the-art explosion simulations from both the sub-M Ch and M Ch channels show that 58 Ni and 56 Ni are not produced in geometric isolation. The forbidden emission lines of Fe II and Ni II in nebular spectra of normal SNe Ia exhibit similar widths and shifts, pointing towards a shared emission region. If indeed 58 Ni and 56 Ni share the volume and excitation conditions then the derived mass ratio of Fe II and Ni II should be representative for Fe/Ni produced in the explosion.
The observed spectra are fit well with our emission model. By using the relation from Section 3.3 we can compute the Ni/Fe ratio. At early times the ratio is still evolving with time as not all the 56 Co has decayed to 56 Fe yet. At late times (>250 d) the ratio remains constant. We find a large spread of Ni/Fe ratios, ranging from 0.02 to 0.08 within the 95 per cent confidence interval. We do not find Figure 7. Comparison between the Fe II and Ni II dominated regions in the optical at 7200 Å (top panel) and the NIR at 20 000 Å (bottom panel) for four observations of SN 2015F and one spectrum of SN 2017bzc. The Ni II lines at 7378 and 19 390 Å originate in the same upper level a 2 F 7/2 and have therefore a fixed line strength ratio that only depends on the ratio of their transition rates. For the atomic data adopted in this work the ratio of the 19 390 Å to the 7378 Å line is 0.202. In the plots, the ratio of the two Ni II emission features is different because of three effects: The optical Ni II feature is a blend of several lines, the flux density is lower at longer wavelengths, and the ratio depends on Galactic as well as host galaxy reddening. Regions of extremely low atmospheric transmission are shaded in grey.  (Seitenzahl et al. 2013, green), the sub-M Ch CO detonation model 'det 2010 1.06' (Sim et al. 2010, orange), and the M Ch 'W7 Z ' model (Nomoto & Leung 2018, red). The black line is not a fit to the data and represents the M Co /M Fe ratio assuming only radioactive decay from 56 Ni to 56 Co to 56 Fe. Grey data points are from Flörs et al. (2018). Black data points are from the newly published objects in this work (SN 2015F andSN 2017bzc). The bottom panel shows the residuals normalized to the pure 56 Ni to 56 Co to 56 Fe decay. Figure 9. Inferred mass ratio of Ni and Fe from archival optical spectra and XShooter observations. Orange data points indicate that the Fe II NIR/VIS ratio was directly inferred from a spectrum covering 4000-20 000 Å. Blue data points indicate optical nebular phase spectra that have been modelled using the relation from Fig. 3 as a prior. Errorbars only indicate the statistical uncertainty from the fit. We assume a rise time of ∼18 d (Ganeshalingam et al. 2011) to compute the time after explosion. The shaded bands display predictions of the Ni to Fe mass ratio from explosion model simulations (Seitenzahl et al. 2013;Shen et al. 2018). Inferred and predicted mass ratios were scaled to t → ∞.
any objects for which we can exclude the contribution of Ni to the nebular phase spectrum.
Our results are in good agreement with sub-M Ch explosions of solar-to supersolar metallicity progenitors. Only few objects have a Ni to Fe ratio that is consistent with explosion predictions from zero-metallicity sub-M Ch WDs. There are only few calculations of non-zero metallicity sub-M Ch explosions Shen et al. 2018). Our data are consistent with both sub-M Ch detonations and double detonations, but they do not allow us to distinguish between these two scenarios.
We find a few objects which have Ni/Fe abundances consistent with nucleosynthetic predictions of exploding M Ch WDs. However, we do not find separate populations but instead the distribution displays a tail of objects which have high Ni/Fe abundances. The abundance distribution of objects which have nebular phase observations peaks at M Ni /M Fe = 0.034 with an 68 per cent Figure 10. The distribution of the Ni/Fe ratio at t → ∞ from all available nebular phase spectra (see Table B1). The Ni/Fe ratio from only optical spectra was corrected according to Section 4.2.1 by σ sys = −0.0033 +0.0037 −0.0041 . The results for SNe with multiple observations were combined so that every supernova in the sample contributes equally to the shown distributionirrespective of the number of spectra. For each unique SN we drew 100 000 samples from the posterior distribution of the Ni/Fe mass ratio. The orange shaded region indicates the region containing 68 per cent of the posterior probability density. The shaded bands display predictions of the Ni to Fe mass ratio from sub-M Ch (Shen et al. 2018, left) and M Ch (Seitenzahl et al. 2013, right)  confidence region between 0.015 and 0.053. Eighty-five per cent of the total probability density falls within the shaded band of sub-M Ch explosion predictions. Our resulting distribution of the Ni/Fe abundance agrees well with the results of Kirby et al. (2019), who determined the Ni/Fe abundance from stellar populations of dwarf galaxies. Only 11 per cent of the total probability lies in the range of M Ch delayed-detonation predictions. The presence of both channels is in agreement with findings from nearby SN remnants (Seitenzahl et al. 2019).
For sub-M Ch we can compare our resulting distribution to explosion yields of progenitors with different masses and metallicities. Progenitors with masses of 0.9 M or less do not produce enough 56 Ni (< 0.3 M ) to explain the brightness of normal SN Ia and are thus discarded for this comparison. The overlap between the range of yields from 1.0 to 1.1 M progenitors with our inferred Ni/Fe distribution is shown in Fig. 10. We find good agreement with progenitors between 0.5 and 2 Z .

C O N C L U S I O N S
The 7200 Å feature in nebular spectra of SNe Ia is composed of emission from Fe II and Ni II and is present in all objects for which this wavelength region has been observed. The relative contributions of the two ions to the feature vary between different SNe. We have presented a method that allows us to place prior constraints on the N e and T and applied it to more than 100 optical archival spectra allowing us to determine the distribution of the Ni/Fe ratio for all objects in our sample. Our main results are as follows: (i) The Fe II emission in the nebular phase can be described by purely thermal forbidden line emission, and it is in agreement with an expanding and cooling nebula.
(ii) The strongest [Fe II] lines in the NIR and at optical wavelengths evolve with time, and the evolution seems to be very homogeneous across our sample. We obtained a relation that describes the evolution of this line ratio. The ratio does not depend on the atomic data. The evolution of the Fe II lines can be used to test more sophisticated spectral synthesis calculations of explosion model predictions -spectra that have been computed from explosion models need to be able to reproduce this relation.
(iii) The 7200 Å feature only contains Fe II and Ni II in normal SNe Ia. A contribution of [Ca II] to this feature would have to be very limited in strength. We used the 19 390 Å line to constrain the 7378 Å line for SN 2015F and SN 2017bzc as these two lines originate from the same upper level. We find no evidence that Ca II] emission is required to reproduce the 7200 Å feature.
(iv) For all objects in the extended sample of more than 100 nebular phase spectra we find that the lines of singly ionized Fe II and Ni II have similar widths and shifts and thus come from the same emitting region and share the same physical conditions. For objects for which NIR spectra are available we can extend this claim to Co II as well.
(v) The display of 130 nebular phase spectra shows a large variety in the relative strengths of the Fe II and Ni II lines in the 7200 Å feature. Translating the relative line strengths into a mass ratio of the singly ionized species results in a distribution which is expected from mainly sub-M Ch explosions.
(vi) We do not find separate populations of sub-M Ch and M Ch explosions. However, the high abundance tail of the distribution extends into the M Ch regime. Eleven per cent of the total probability distribution lies within the M Ch predictions of the Ni/Fe abundance. The spectra have been corrected for redshift and Galactic extinction to better illustrate the position of the strongest iron, nickel, and cobalt lines (dashed vertical lines). Additional information on these observations can be found in Table B1.

A P P E N D I X B : OV E RV I E W O F N E B U L A R S P E C T R A
In Table B1, we provide the SN name, subtype, combined Galactic and host galaxy colour excess, redshift, and the date of B-band maximum for each SN Ia that is used in the analysis. Multiple observations of the same SN Ia are sorted by increasing epoch. We also show the telescope and instrument that was used to obtain the spectrum. The measured Ni/Fe mass ratio in the limit t → ∞ is also given in Table B1 for each spectrum.