Including a Luminous Central Remnant in Radiative Transfer Simulations for Type Iax Supernovae

Type Iax supernovae (SNe Iax) are proposed to arise from deflagrations of Chandrasekhar mass white dwarfs (WDs). Previous deflagration simulations have achieved good agreement with the light curves and spectra of intermediate-luminosity and bright SNe Iax. However, the model light curves decline too quickly after peak, particularly in red optical and near-infrared (NIR) bands. Deflagration models with a variety of ignition configurations do not fully unbind the WD, leaving a remnant polluted with $^{56}\mathrm{Ni}$. Emission from such a remnant may contribute to the luminosity of SNe Iax. Here we investigate the impact of adding a central energy source, assuming instantaneous powering by $^{56}\mathrm{Ni}$ decay in the remnant, in radiative transfer calculations of deflagration models. Including the remnant contribution improves agreement with the light curves of SNe Iax, particularly due to the slower post-maximum decline of the models. Spectroscopic agreement is also improved, with intermediate-luminosity and faint models showing greatest improvement. We adopt the full remnant $^{56}\mathrm{Ni}$ mass predicted for bright models, but good agreement with intermediate-luminosity and faint SNe Iax is only possible for remnant $^{56}\mathrm{Ni}$ masses significantly lower than those predicted. This may indicate that some of the $^{56}\mathrm{Ni}$ decay energy in the remnant does not contribute to the radiative luminosity but instead drives mass ejection, or that escape of energy from the remnant is significantly delayed. Future work should investigate the structure of remnants predicted by deflagration models and the potential roles of winds and delayed energy escape, as well as extend radiative transfer simulations to late times.


INTRODUCTION
Through modern transient surveys it has become clear that Type Ia supernovae (SNe Ia) are a diverse class of objects: in addition to the "normal" SNe Ia, which exhibit an empirical relation between their peak brightness and light curve evolution (Phillips 1993), there are multiple distinct sub-classes (see e.g.Taubenberger 2017 for a review).
Type Iax supernovae (SNe Iax, Li et al. 2003;Foley et al. 2013) are the most numerous peculiar sub-class of SNe Ia, recently estimated by Srivastav et al. (2022) to make up 15 +15 −9 % of the total SNe Ia rate.Compared to normal SNe Ia, SNe Iax show low absolute brightness at peak relative to their light curve widths (e.g.McClelland et al. 2010;Stritzinger et al. 2014) and thus fall outside the widthluminosity relation.They display a wide diversity in their band light curve properties showing variations in peak brightness of up to six magnitudes in certain bands and diversity in their light curve rise times and decline rates (Magee et al. 2016;Srivastav et al. 2022).
SNe Iax are spectroscopically distinct from normal SNe Ia: the ★ E-mail: f.callan@qub.ac.uk blue continua and absorption from higher ionisation species observed in their early-time spectra mean they show more similarities to the sub-classes of 1991T-like and 1999a-like SNe Ia, which exhibit hot photospheres (Foley et al. 2013;McCully et al. 2014).Additionally, SNe Iax spectra show significantly lower expansion velocities at peak compared to normal SNe Ia (Foley et al. 2013;Stritzinger et al. 2015).
The late time spectra of SNe Iax show a large number of permitted lines from iron group elements (IGEs) along with some intermediate mass elements (IMEs) and potential signatures of oxygen, while forbidden lines of IGEs and IMEs are also present (Jha et al. 2006;Foley et al. 2016).As such, SNe Iax never enter a fully nebular phase (unlike normal SNe Ia Graham et al. 2017;Taubenberger et al. 2015).
The SNe Iax sub-class also shows a large spectroscopic diversity in line velocities, strengths and widths particularly at later times (McCully et al. 2014;Stritzinger et al. 2015;Yamanaka et al. 2015;Foley et al. 2016).
Multiple explosion scenarios have been suggested to explain SNe Iax, such as the core collapse of massive stars (Foley et al. 2009;Valenti et al. 2009;Moriya et al. 2010), pulsationally delayed detonations (Stritzinger et al. 2015), the failed detonation of an oxygen-neon (ONe) WD merging with a carbon-oxygen (CO) secondary (Kashyap et al. 2018), a WD-neutron star merger (Fernández & Metzger 2013;Bobrick et al. 2022) and a detonation in a helium shell surrounding a CO WD (Foley et al. 2009).
However, amongst the most widely discussed scenarios is that in which SNe Iax are pure deflagrations of Chandrasekhar mass ( Ch ) CO WDs (Branch et al. 2004;Phillips et al. 2007), as the properties of such explosions naturally reproduce many of the observed characteristics of SNe Iax.In particular, simulations of pure deflagrations in  Ch CO WDs (Jordan et al. 2012;Kromer et al. 2013Kromer et al. , 2015;;Fink et al. 2014;Long et al. 2014;Leung & Nomoto 2020;Lach et al. 2022) reproduce many of the observed characteristics of SNe Iax, such as lower 56 Ni masses and hence sub-luminous light curves as well as lower expansion velocities.Furthermore, the pure deflagration model sequences presented by Kromer et al. (2013); Fink et al. (2014) and Lach et al. (2022) (hereafter K13, F14 and L22 respectively) are together able to reproduce almost the full diversity in peak brightness of observed SNe Iax by varying the initial geometry of the explosion simulations, failing to match only the faintest observed SNe Iax.These models also produce light curves and spectra in reasonably good agreement with bright and intermediate-luminosity SNe Iax, particularly at earlier epochs.Additionally, Camacho-Neves et al. (2023) have demonstrated good spectroscopic agreement between a deflagration model and the SN Iax, SN 2014dt, over ∼ 500 days of its evolution.
The near-infrared light curves of SNe Iax further support the pure deflagration scenario: the lack of a prominent secondary maximum that is observed for normal SNe Ia (Li et al. 2003) points to a wellmixed ejecta structure for SNe Iax (Jha et al. 2006;Phillips et al. 2007;Sahu et al. 2008), consistent with turbulent deflagration models (although see Stritzinger et al. 2015 andMagee et al. 2022 for further discussion of stratification in SNe Iax).
Despite the generally good agreement between  Ch CO WD pure deflagration models and observed SNe Iax, some systematic differences remain.Specifically, the K13, F14 and L22 model light curves decline too quickly post peak, particularly in the red and NIR bands, meaning the bolometric light curves also decline too rapidly.This systematic difference worsens when comparing to fainter SNe Iax where the decline post peak of the models is too fast compared to observed SNe Iax in all bands.The faster decline for red, optical, and NIR wavelengths after peak relative to other bands also makes the model spectra too blue at later times relative to observed SNe Iax.
The outcome of pure deflagration simulations is sensitive to the initial choice of ignition conditions (Jordan et al. 2012;Kromer et al. 2013;Fink et al. 2014;Long et al. 2014;Kromer et al. 2015;Leung & Nomoto 2020;Lach et al. 2022).For a wide variety of ignition configurations, the explosion is not energetic enough to fully unbind the WD (Jordan et al. 2012;Kromer et al. 2013;Fink et al. 2014) potentially leaving behind a stellar remnant polluted with burning products from the deflagration, including 56 Ni.The remnant is therefore expected to be luminous.Such a luminous remnant may have been detected in late time observations (∼ 4 years after explosion) of the faint SN Iax, SN 2008ha (Foley et al. 2014) and bright SN Iax, SN 2012Z (McCully et al. 2022).Further evidence of a remnant may be present in the slowly declining late time light curves of SN 2014dt (Kawabata et al. 2018) and SN 2019muj (Kawabata et al. 2021) and the peculiar late time spectra of SNe Iax that never become truly nebular (Kromer et al. 2013;Foley et al. 2016;Kawabata et al. 2018Kawabata et al. , 2021;;Maeda & Kawabata 2022;Camacho-Neves et al. 2023).
Accounting for the remnant may be key to fully understanding the observed properties of SNe Iax.An exploratory calculation by K13 found adding the energy deposited by radioactive 56 Ni in the remnant to the light curve predicted for their N5def deflagration model significantly improved the agreement with the bolometric light curve of SN 2005hk, a bright SNe Iax, particularly at later times.However, their calculation simply adds this energy to the lightcurve computed from the ejecta model and does not take into account the interaction of this radiation with the ejecta.Therefore, with this approach only a crude bolometric light curve can be calculated and the early-time bolometric light curves will not be accurately represented.
Shen & Schwab (2017) also found their models of post SNe Iax remnant winds to be consistent with late time photometric observations of SNe Iax.Additionally, from analysis of the late-time spectra of the intermediate-luminosity SNe Iax, SN 2019muj, Maeda & Kawabata (2022) find that a denser inner ejecta component is present.They attribute this to a secondary mass ejection from a remnant.The faint hybrid carbon-oxygen-neon (CONe) WD pure deflagration model which Kromer et al. (2015) produced as an attempt to match SN 2008ha also predicts a luminous remnant which they suggest may have a significant impact on the optical display of the model.
The aim of this paper is to incorporate a central source representing the remnant material into our 3D, time-dependent radiative transfer simulations in order to predict how a luminous remnant impacts the light curves and spectra predicted for  Ch CO WD pure deflagration models covering a range of peak luminosities.This improves on the exploratory calculation of K13 since we can simulate the early time diffusion of this radiation through the ejecta and investigate how it affects light curves, colours and spectra.We take the pure deflagration models presented by L22 as a starting point.Informed by the remnant compositions predicted by the hydrodynamic explosion simulations of L22 we add a central 56 Ni source representing a luminous remnant to a sub-set of the standard L22 deflagration models.We discuss the treatment we adopt for the remnant and summarise our radiative transfer simulation set up in Section 2. In Section 3 we make detailed comparisons between the light curves and spectra predicted by our radiative transfer simulations and those of observed SNe Iax, before discussing the implications of our results along with our conclusions in Section 4.

Explosion models
The sequence of  Ch CO WD pure deflagration models presented by L22 are the basis for this work.These models all adopt a singlespark ignition, which has been argued to be a more realistic ignition configuration than a multi-spark ignition as adopted by F14 (Kuhlen et al. 2006;Zingale et al. 2009Zingale et al. , 2011;;Nonaka et al. 2012).The L22 sequence provides a good starting point to explore the impact of including a luminous remnant in our simulations: all models in this sequence produce explosions that are not energetic enough to fully unbind the WD, leaving behind remnants with masses greater than a solar mass that are polluted with burning products from the deflagration including 56 Ni.Additionally, the L22 model sequence covers almost the full peak brightness range of the SNe Iax class, failing to reproduce only the very brightest and faintest observed SNe Iax.We can therefore explore how a luminous remnant impacts our model comparisons with a wide variety of observed SNe Iax from bright to faint events.
The leafs (Reinecke et al. 1999) hydrodynamic explosion simulations used to calculate the L22 models adopted an expanding grid to track the ejecta (Röpke 2005).This means the remnant is not well resolved at the end of the hydrodynamic explosion simulations (at  = 100 s).Therefore, only basic properties of the remnant can be extracted: the overall mass and composition at the end of the explosion simulations can be determined from the tracer particles which remain bound (see L22 for details).We note however, that while we can say the remnant is bound at  = 100 s, this does not necessarily mean that the material in the remnant will all remain bound at later times.Indeed, the energy predicted to be deposited in the remnants from the 56 Ni decay chain indicates they are super-Eddington by a factor of ∼ 1000 to 10000 depending on the model and time in the simulation.This means it is highly likely that some or all of the remnant material becomes unbound, perhaps in the form of a radioactively driven post supernova wind as proposed by Shen & Schwab (2017) or as a secondary mass ejection as described by Maeda & Kawabata (2022).Therefore, when we discuss the remnant in this work we are referring to any material which remains bound up to the end of the hydrodynamic explosion simulation at 100 s and so would not be included as part of the ejecta in the radiative transfer simulations for the standard L22 model sequence.
The remnants' structure at the end of the L22 explosion simulations is that of a puffed-up envelope of burning products from the deflagration, mixed with CO material, which is settling onto a relatively dense CO core (see also K13 and Jordan et al. 2012).The remnant contains a significant fraction of 56 Ni (between 2 and 7% depending on the model) with this 56 Ni predicted to be primarily found in the envelope (again see L22 for details).For the models in this work, we calculated an estimate of the maximum velocity the remnant material would approach if all the radioactive energy deposited in the remnant over the simulation time was converted into kinetic energy shared between all the remnant mass.The estimated velocities (between ∼ 500 to 800 km s −1 ) are comparable to the velocity of a single grid cell for all models.However, it is likely that only the 56 Ni rich part of the remnant envelope is ejected.The Shen & Schwab (2017) remnant model with most similar mass to our model remnants has an envelope mass of 0.1 M ⊙ (see Section 4.2 for more detailed discussion).If we assume all the radioactive energy is deposited in the remnant envelope during the simulation time, the envelope would reach an estimated velocity between ∼ 2000 to 2600 km s −1 depending on the model.While more significant, these velocities are still only comparable to the velocities of a small number of inner grid cells in our simulations.As such, even if our model remnants become partially or fully unbound they will likely remain small on the scale of the computational grid.

Incorporating the remnant in artis
To calculate synthetic observables (light curves and spectra) for our models, to compare with observed SNe Iax, we use the 3D timedependent Monte Carlo radiative transfer code artis 1 (Sim 2007;Kromer & Sim 2009;Bulla et al. 2015;Shingles et al. 2020).In the standard mode of operation, artis follows the propagation of -ray photons emitted by radioactive decays and deposits energy in the supernova ejecta before solving the radiative transfer problem selfconsistently.For this study we also include the contribution from a luminous remnant source in our radiative transfer simulations by placing a 56 Ni source in the centre of the models with 56 Ni mass informed by the remnant 56 Ni masses predicted by the L22 hydrodynamic explosion simulations.The main properties of the remnant and ejected material predicted by the L22 explosion simulations that form the basis of our new models are summarised in the 56 Ni decay chain in the remnant contributes to the optical display of the models the remnant will have a larger impact on the observed properties of faint models.artis is not configured to deal with a central luminous source in its standard mode of operation.Therefore, to incorporate a treatment for luminous remnants, we made use of developments made by Collins et al. (2023a) which allows more flexibility in the way energy is injected in our artis simulations.We extended this method to treat energy due to radioactive decays from the 56 Ni decay chain in the remnant separately to energy injected by radioactive decays in the ejecta.The rate and amount of radiation escaping the remnant is affected by a number of factors including the diffusion of the radiation through the remnant and potential mechanisms driven by the radioactive energy deposited in the remnant such as winds in the remnant envelope (Shen & Schwab 2017) or the ejection of mass from the remnant (Maeda & Kawabata 2022).However, for simplicity, in our prescription we assume the remnant emits instantaneously according to the 56 Ni decay chain (i.e.assuming no significant delay due to diffusion in the remnant itself).Due to the compactness of the remnant material we expect the gamma rays emitted by the remnant to be trapped with their energy used to heat the remnant (see also discussion in K13).We therefore assume the remnant emits as a black body.Specifically, each time there is a 56 Ni decay in the remnant, a packet of UVOIR radiation is emitted from the centre of the model grid and assigned a co-moving frame frequency sampled from a black body distribution with the temperature treated as a parameter of the simulation and representing the remnant temperature at the decay time (see discussion below).Once emitted, the UVOIR radiation packets travel through the ejecta and interact in the same way as UVOIR radiation packets originating in the ejecta.This approach allows us to obtain band light curves and spectra for our pure deflagration models with the contribution of a luminous remnant included, which we can compare directly to the light curves and spectra of observed SNe Iax.
In this work, we explore simulations first assuming a constant remnant temperature.We also test the sensitivity of our findings to the remnant temperature treatment adopted by exploring a temperature evolution of the remnant consistent with black body emission from a surface of constant radius.The SNe Iax we compare to throughout this study have temperature estimates from observations which provide an approximate lower limit on the temperatures of potential remnants present for these objects.These temperature estimates guide the range of remnant temperatures we explore in this work: Sahu et al. (2008)  our models have estimated temperatures of ∼ 3000 K to 6000 K at 100 days after explosion.Throughout this paper, the standard pure deflagration models retain their names as introduced by L22 with the structure rX_dY_Z where r refers to the ignition spark offset from the centre of the WD in km, d is the WD central density in 10 9 g cm −3 and Z the metallicity relative to the solar value (all models presented in this work adopt solar metallicity).We differentiate the models including a remnant contribution with an "R" in the model name and also include information on the remnant temperature: for models with fixed remnant temperature we include this temperature in Kelvin in the model name and differentiate our model which adopts a temperature evolution consistent with a black body of constant radius with the "const_rad" label.Models which adopt remnant 56 Ni masses less than those predicted by the L22 hydrodynamic explosion simulations of the models are labelled with the fraction of the predicted 56 Ni mass that is adopted for the remnant.For example model r48_d5.0_Z_R_6000K_0.33Ni is based on the standard r48_d5.0_Zpure deflagration model but includes the contribution from a remnant with a fixed characteristic temperature of 6000 K and a 56 Ni mass a third of that predicted by the L22 hydrodynamic explosion simulation of the standard deflagration model.We summarise the assumptions of the remnants for the models presented in this work in Table 2.

Radiative transfer
artis follows the methods of Lucy (2002Lucy ( , 2003Lucy ( , 2005) ) and is based on dividing the radiation field into indivisible energy packet Monte Carlo quanta.In this work we adopt the NLTE (non local thermodynamic equilibrium) approximation described by Kromer & Sim (2009) and as used by F14 and L22.We do not make use of the full NLTE and non-thermal capabilities added to artis by Shingles et al. (2020) which are required to accurately model the nebular phase of SNe Ia (see, e.g., Dessart et al. 2014;Jerkstrand 2017;Shingles et al. 2022).As such we do not extend the simulations presented in this work to nebular phases and halt our analysis at 80 days post explosion.We mapped the 3D abundance and density structure of the unbound material at the end of the L22 hydrodynamic simulations (by which point homologous expansion is a good approximation Röpke 2005) to a 50 3 Cartesian grid.3D radiative transfer simulations were then carried out for each model.In each simulation, 3×10 7 energy packets were tracked through the ejecta for 200 logarithmically spaced time steps between 0.3 and 100 days.We used the atomic data set compiled by Gall et al. (2012) which is sourced from Kurucz & Bell (1995) and Kurucz (2006).The much more extensive line list of Kurucz (2006) is used for the important second and third ionisation stages of Fe, Co and Ni with the Kurucz & Bell (1995) atomic data used for the remaining ions.We adopted a grey approximation in optically thick cells (c.f.Kromer & Sim 2009) and at times earlier than 0.4 days post explosion, local thermodynamic equilibrium (LTE) was assumed.

RESULTS
In the following section, we present comparisons between the bolometric and band light curves and spectra predicted by radiative transfer simulations of our new models with remnant contribution included and the standard L22 deflagration models along with observed SNe Iax.To investigate the impact of including the remnant over a wide range of model brightnesses we include comparisons between bright, intermediate-luminosity and faint models along with corresponding observed SNe Iax (the light curve properties of our new models and the L22 models we compare to are summarised in Tables A1 and A2).The intermediate-luminosity and faint models with remnant contribution included that produced best agreement with observed SNe Iax adopt remnant 56 Ni masses less than those predicted by the hydrodynamic explosion simulations of the models.We discuss this in detail in Section 4. L22 found noticeable viewing angle effects for their standard pure deflagration models but these effects were not so large that they altered any of their findings.This is consistent with the viewing angle effects usually exhibited by pure deflagration models which are generally less substantial than those of other SNe Ia models due to their well-mixed ejecta structures (see e.g.K13, F14 and L22).We find our models including central remnants do not show significant viewing angle effects other than those already exhibited by the standard L22 pure deflagration models.Therefore, we focus on presenting angle-averaged results since these illustrate our main findings, but we comment on observer orientation effects where relevant.We also include the range in model parameters obtained from the line-of-site dependent light curves alongside the angle averaged values in Tables A1 and A2 for reference.

Light curves
Figure 1 shows comparisons between the bolometric light curves of the brightest standard pure deflagration model from the L22 sequence, our new r10_d4.0_Z_R_8000Kmodel and the bright SN Iax, SN 2005hk.Also plotted is the bolometric light curve produced if energy deposited by radioactive decays in the remnant is instantaneously radiated (labelled "Energy_dep").Following the approach of Kromer et al. (2013) we also plot the bolometric light curve produced if we simply add this remnant energy deposition to the bolometric light curve of the r10_d4.0_Zmodel (labelled r10_d4.0_Z+ Energy_dep).From Figure 1 we see the r10_d4.0_Z+ Energy_dep light curve produces good agreement with the bolometric light curve of SN 2005hk at later times (after ∼ 50 days).However, at early times it significantly overestimates the brightness of the bolometric light curve until ∼ 10 days after explosion and then declines too quickly until ∼ 40 days.The r10_d4.0_Z_R_8000Kbolometric light curve provides a good match to both the rise, decline and peak bolometric brightness of the SN 2005hk light curve, showing significantly improved agreement compared to the r10_d4.0_Z+ R_energy_dep bolometric light curve at earlier times.Therefore, including the energy injected in the remnant directly in the radiative transfer simulations is important for accurate predictions of the remnant contribution to the early time bolometric light curves.Relative to Model r10_d4.0_Z the r10_d4.0_Z_R_8000Kbolometric light curve provides a better match to SN 2005hk in rise and particularly in decline after peak also producing much better agreement with the peak brightness of SN 2005hk.We also note that the brightest viewing angles of Model r10_d4.0_Z_R_8000Kreach the peak bolometric luminosity of the brightest members of the SNe Iax class within their uncertainties, something not previously achieved by the L22 pure deflagration models that adopt a single-spark ignition.
Figure 2 shows angle-averaged bolometric and UBVRĲH-band light curves for the r10_d4.0_Zmodel and four new models from this work in which energy injected by radioactive decays in the remnant are also included in the radiative transfer.The light curves of SN 2005hk are included for comparison.Each of the four models with remnant contribution included adopt the full remnant 56 Ni mass predicted by the hydrodynamic explosion simulations of L22.However, the remnant temperature treatment is different for each model so we can investigate the sensitivity of the synthetic observables obtained from our radiative transfer simulations to what is assumed about the remnant's emission.The r10_d4.0_Z_R_2000K,r10_d4.0_Z_R_8000Kand r10_d4.0_Z_R_15000Kmodels all have fixed characteristic black body temperatures while the r10_d4.0_Z_R_const_radmodel corresponds to the remnant emitting as a black body surface of fixed radius (i.e. the remnant temperature decreases with time).We discuss how these different choices of remnant temperature treatment impact the band light curves of the models in detail below.However, we note that provided the energy is injected centrally in the remnant, the optical light curves are insensitive to different choices of remnant SED until ∼ 30 days after explosion (see Figure 2).
From Figure 2 we see that Model r10_d4.0_Z is not bright enough at peak and also too fast in both rise and decline in all bands to match the light curves of SN 2005hk.By comparison the r10_d4.0_Z_R_8000K,r10_d4.0_Z_R_15000Kand r10_d4.0_Z_R_const_radmodels provide a significantly better match to the peak brightness, rise and decline of the band light curves of SN 2005hk (as was the case for the bolometric comparisons).Specifically, we emphasise the extremely good agreement between the R, I and H-band model light curves and those of SN 2005hk, in particular in decline after peak.Compared to bright SNe Iax such as SN 2005hk, previous pure deflagration models (e.g. the N5def model presented by K13) declined significantly too quickly after peak in these bands even when the decline was well matched in blue bands.Including the remnant contribution in our simulations significantly improves the agreement with the red optical and NIR bands although there is still a deficit in J-band albeit much reduced.
The band light curves of the r10_d4.0_Z_R_8000Kand r10_d4.0_Z_R_15000Kmodels are very similar (see Figure 2).The biggest differences between the two models is in J and H-band where the r10_d4.0_Z_R_8000Kmodel produces better agreement with SN 2005hk.This marginal improvement in agreement does not, however, rule out the higher remnant temperature of 15000 K. Indeed, through additional simulations we have confirmed that a wide range of fixed remnant temperatures from 6000 K to 15000 K (the highest temperature we tested) produce model light curves in good agreement with those of SN 2005hk.Therefore, our results are not excessively impacted by different choices of fixed remnant temperature.
Model r10_d4.0_Z_R_const_radallows us to test the impact of instead assuming the remnant radius is fixed, meaning the remnant temperature evolves in time.For Model r10_d4.0_Z_R_const_radwe adopt a remnant radius of 4.9 × 10 14 cm.This is broadly consistent with the photospheric radius Maeda & Kawabata (2022) derive for the proposed inner ejecta component of SN 2019muj, an intermediate luminosity SN Iax.The remnant radius of Model r10_d4.0_Z_R_const_radyields a black body temperature of 8000 K at 30 days after explosion (i.e. the time the optical bands become sensitive to the remnant temperature).We note that formally this radius is too large to be consistent with central energy injection in the remnant for early times.Nevertheless we continue to assume central energy injection noting the emerging optical light curves are insensitive to the assumed remnant SED prior to this time.Indeed, by the time the choice of remnant SED becomes relevant, it is already contained within the inner 15% of the ejecta of Model r10_d4.0_Z_R_const_radwhich is consistent with our assumption of central energy injection in the remnant.This is true for all our models presented in this study that produce a reasonable match to observed SNe Iax.
Comparing band light curves for the r10_d4.0_Z_R_const_radand r10_d4.0_Z_R_8000Kmodels, we see adopting an evolving remnant temperature instead of a fixed temperature leads to only very small variations in the band light curves, even at later times.The r10_d4.0_Z_R_const_radmodel shows a remnant effective temperature variation from ∼ 12500 K to 6000 K over the simulation time from 0.3 to 100 days.Our results are therefore not overly sensitive to remnant temperature variations in this range over the simulation time.The most noticeable differences are observed in J and H-bands where the r10_d4.0_Z_R_const_radmodel is slightly brighter than the r10_d4.0_Z_R_8000Kmodel from ∼ 50 days onwards.This is because the r10_d4.0_Z_R_const_radmodel has lower effective remnant temperature at later times which increases the average wavelength of photon packets emitted by the remnant and thus the brightness in NIR bands.Overall, any variations in the band light curves between the r10_d4.0_Z_R_8000Kand r10_d4.0_Z_R_const_radmodels are not significant in relation to comparisons made with the light curves of SN 2005hk even in J and H-bands.Therefore, within the time range of our simulations our conclusions are not sensitive to whether the remnant has a fixed or evolving temperature and we make the choice to primarily use remnants with fixed temperature for this study.
Of the models with remnant contribution included Model r10_d4.0_Z_R_2000Kproduces by far the poorest agreement with SN 2005hk, declining much too quickly post peak in the optical (2014) of ≲ 4500 K and ∼ 4500 to 7000 K respectively.These estimates do not, however, rule out the remnant temperatures adopted for the remaining models we compare to SN 2005hk in Figure 2.

Spectra
Figure 3 shows spectral comparisons, at phases relative to Bpeak, in absolute flux, between Model r10_d4.0_Z_R_8000K, the r10_d4.0_ZL22 model and SN 2005hk.Also plotted are the r10_d4.0_Zspectra scaled by a constant factor across all wavelengths such that the peak of the r10_d4.0_Zspectra matches the peak of the r10_d4.0_Z_R_8000Kmodel spectra at each epoch compared.All spectra scaled throughout this work follow this scaling approach.We compare relative to B-peak so the flux of the models is as similar as possible across the epochs compared.Compared to the r10_d4.0_Zmodel, the r10_d4.0_Z_R_8000Kmodel spectra provide a significantly better match to the absolute flux of the SN 2005hk spectra across all epochs.This was already clear from the band light curve comparisons discussed above, but there are also some interesting differences in spectral shape and features, which we discuss below.The r10_d4.0_Z_R_8000Kmodel shows its poorest agreement with the spectra of SN 2005hk when comparing at the earliest epoch 8.4 days before B-peak.Although the SED of the r10_d4.0_Z_R_8000Kmodel provides a reasonable match to SN 2005hk, particularly for wavelengths ≳ 5500 Å, the spectrum is too faint.This is a result of the rapid early-time light curve evolution.When the spectra are instead compared at epochs relative to explosion (corresponding to a spectral epoch 1.3 days later than that plotted for Model r10_d4.0_Z_R_8000K) the flux match is significantly improved for this earliest epoch.From Figure 3 we also see the r10_d4.0_Z_R_8000Kmodel spectrum shows a more noticeable lack of flux for wavelengths ≲ 4000 Å, with the r10_d4.0_Zscaled model producing better agreement with SN 2005hk over these wavelengths.This is a result of the shorter rise time to B-peak of the r10_d4.0_Zmodel.This means the r10_d4.0_Zmodel is closer to explosion for the comparison shown and thus has hotter ejecta and therefore a bluer spectrum.When the models are instead compared relative to explosion they show extremely similar spectral evolution until ∼ 10 days since explosion.For this earliest epoch both models also show stronger spectral features for wavelengths ≲ 5500 Å when comparing relative to explosion, improving the spectral agreement with SN 2005hk.
For the epochs at 4.3 days before B-peak and 3.6 and 14.5 days after B-peak the r10_d4.0_Z_R_8000Kand r10_d4.0_Zscaled models display spectra which are very similar across all wavelengths.Both models display a reasonable match to the SED of SN 2005hk across these epochs.For the epochs at 4.3 days before and 3.6 days after B-peak the agreement relative to the earliest epoch is primarily improved for wavelengths ≲ 5500 Å.However, at 14.5 days after B-peak, the epoch the models show best overall spectroscopic agreement with SN 2005hk, there is an improved match to the SED of SN 2005hk across all wavelengths.Additionally, over these epochs, both models show improved agreement with the spectra of SN 2005hk relative to the earliest epoch shown due to the increased number and strength of spectral features.At 23.7 days after B-peak, the model spectra remain very similar although the r10_d4.0_Z_R_8000Kmodel matches the strengths of some of the absorption features in the blue wavelengths observed for SN 2005hk better, such as the strong absorption trough at ∼ 5200 Å primarily due to Fe ii.While the overall spectral agreement with SN 2005hk is still good for this epoch, both models show emission features (primarily attributed to Fe and Co ii in the models) that are too weak relative to SN 2005hk, particularly for red wavelengths.For the earlier epochs shown, the model spectral features are in general more blue-shifted than those of SN 2005hk.However, for later epochs, although there are still some mismatches, the model spectral features do not show a preference for being blue-shifted or red-shifted relative to SN 2005hk.This may suggest the outer regions of the model ejecta, where spectra are formed at early times, have velocities in general too high compared to SN 2005hk, while the inner ejecta of the models, have velocities which provide a better match to SN 2005hk.Another explanation for this behaviour may be differences in the ionisation state between the outer ejecta of SN 2005hk and the models leading to the early-time spectral features of the models generally being formed at higher velocities.
The largest difference between the r10_d4.0_Z_R_8000Kand r10_d4.0_Zscaled model spectra is found for the latest epoch shown (39.5 days after B-peak).At this epoch, the r10_d4.0_Z_R_8000Kspectrum matches the SED of SN 2005hk well across all wavelengths shown and reproduces the strength and velocities of many of the spectral features across a wide range of wavelengths.The r10_d4.0_Zscaled spectrum, on the other hand, has a noticeable lack of flux for wavelengths longer than ∼ 5500 Å and has spectral features much weaker than SN 2005hk for these redder wavelengths.This lack of flux for red wavelengths at later times is a general issue for standard pure deflagration models (see Section 1).Overall, the r10_d4.0_Z_R_8000Kand r10_d4.0_Zscaled model spectra are very similar for all but the latest epoch shown, producing similarly good spectroscopic agreement with SN 2005hk over these epochs.

Light curves
Figure 4 shows band light curves plotted for Model r48_d5.0_Z_R_6000K_0.33Ni(which adopts a remnant 56 Ni mass a third of that predicted by L22, see Table 2 for summary of remnant properties), Model r114_d6.0_Z_R_4000K, the L22 r48_d5.0_Zmodel and the intermediate-luminosity SN Iax, SN 2019muj.Model r48_d5.0_Z_R_6000K_0.33Ni,our model with remnant contribution included in best agreement with SN 2019muj, produces significantly better agreement with the light curves of SN 2019muj than the closest matching L22 model r48_d5.0_Z,matching the peak, rise and particularly decline of SN 2019muj better in all bands (apart from band, where the small number of photometric points for SN 2019muj makes the level of agreement hard to judge between models).In particular, the remnant contribution removes the systematic problem of the light curves declining too rapidly.The improvement is again most dramatic in the red optical and NIR bands.While the r48_d5.0_Z_R_6000K_0.33Nimodel produces best agreement with the optical band light curves of SN 2019muj, the agreement is still reasonably good for all NIR bands shown, despite the rise being a bit slow relative to SN 2019muj in the H and K s -bands.The remnant temperature adopted for Model r48_d5.0_Z_R_6000K_0.33Ni is consistent with the photospheric temperature of ∼ 5500 K estimated by Maeda & Kawabata (2022) for SN 2019muj at 131 days after explosion.Overall, including energy injected by radioactive decays in the remnant significantly improves the agreement of pure deflagration models with the band light curves of intermediate-luminosity SNe Iax such SN 2019muj.
To match the optical luminosity of SN 2019muj for models including the entire remnant 56 Ni mass predicted by L22 we require remnant temperatures lower than that adopted for Model r48_d5.0_Z_R_6000K_0.33Ni.Of the models we investigated that adopted the entire L22 56 Ni mass, Model r114_d6.0_Z_R_4000Kproduces best agreement with SN 2019muj.Comparing the light curves of this model and Model r48_d5.0_Z_R_6000K_0.33Nito SN 2019muj we see both models produce good agreement with the BVgr band light curves of SN 2019muj although Model r48_d5.0_Z_6000K_0.33Niprovides an overall better match in these bands.Model r48_d5.0_Z_6000K_0.33Nialso provides a reasonably good match to SN 2019muj in all other bands shown.However, Model r114_d6.0_Z_R_4000K is too bright at peak in all remaining bands relative to SN 2019muj.Additionally, all the NIR light curves of Model r114_d6.0_Z_R_4000K are at least a magnitude too bright compared to SN 2019muj from peak until the end of the simulation at 100 days.Model r114_d6.0_Z is the faintest from the L22 sequence meaning it is not possible to reproduce the correct brightness for any observed SNe Iax fainter than SN 2019muj if we assume the entire remnant 56 Ni mass contributes directly to the observed luminosity of the model.Therefore, to produce good agreement with both intermediate-luminosity and faint SNe Iax, we must adopt remnant 56 Ni masses in our models less than those predicted by L22.We discuss the implications of this finding in more detail in Section 4.

Spectra
Figure 5 shows absolute flux spectral comparisons between Model r48_d5.0_Z_6000K_0.33Ni,Model r48_d5.0_Z,Model r48_d5.0_Zscaled to match the peak of r48_d5.0_Z_6000K_0.33Niand SN 2019muj.Model r48_d5.0_Z_6000K_0.33Niprovides a significantly better overall flux match to the spectra of SN 2019muj at all epochs compared to the r48_d5.0_Zmodel.Therefore, following the same approach as Section 3.1.2we primarily focus on the r48_d5.0_Zscaled flux model alongside Model r48_d5.0_Z_6000K_0.33Nifor our comparisons to SN 2019muj so we can evaluate how including the remnant contribution impacts the spectral evolution of intermediateluminosity deflagration models.
The r48_d5.0_Z_6000K_0.33Nimodel provides a very good match to the overall SED and the strengths and velocities of a significant number of the spectral features of SN 2019muj, across all epochs shown.Relative to the r48_d5.0_Z_6000K_0.33Nimodel, the r48_d5.0_Zscaled model produces a similarly good match to the spectral features of SN 2019muj except for the latest epoch shown, where the r48_d5.0_Zspectral features appear too weak for wavelengths longer than ∼ 6500 Å primarily due to the lack of flux for the model at these red wavelengths.The r48_d5.0_Zscaled model produces a significantly worse match to the colour evolution of SN 2019muj: while the model matches the first epoch shown (5 days before B-peak), it becomes too red at 2 and 11 days after peak.At 19.7 days after peak the model then exhibits an improved match to the colours of SN 2019muj before showing spectra which become increasingly too blue for the latest two epochs shown at 27.7 and 50.1 days after peak.Compared to bright models, including the remnant contribution for intermediate-luminosity models leads to a greater change in spectral properties and bigger overall improvement in the agreement with observed SNe Iax.This suggests the energy injected into the ejecta due to the radiation emitted from the remnant has a  (Persson et al. 1998).The K s filter is a K filter which has been modified to reduce thermal background for ground based telescopes.
greater impact on the ejecta conditions of intermediate-luminosity deflagration models.

Light curves
From Figure 4 we can see band light curve comparisons between Model r114_d6.0_Z_6000K_0.1Ni(which adopts a remnant 56 Ni mass a tenth of that predicted by L22), the r114_d6.0_ZL22 model and the faint SNe Iax, SN 2020kyg.From fitting the SED of SN 2020kyg, Srivastav et al. (2022) find a black body temperature which varies from ∼ 8000 K at g band peak to ∼ 4000 K at 60 days after g band peak, which is consistent with the remnant temperature we adopt for our model.L22 predict a remnant 56 Ni mass ∼ 5 times greater than the ejecta 56 Ni mass for Model r114_d6.0_Z.However, the ejecta component of Model r114_d6.0_Z is already brighter than SN 2020kyg at peak in most bands.Therefore, it follows that to provide a reasonable match to the peak brightness of SN 2020kyg our model with remnant contribution included must adopt a 56 Ni mass only a small fraction of that predicted by L22 (see Section 4 for discussion).The r114_d6.0_Z_6000K_0.1Nimodel matches the decline of SN 2020kyg in all bands better than the standard r114_d6.0_Zmodel, which declines much too rapidly in all bands.In particular the r114_d6.0_Z_6000K_0.1Nimodel quite successfully matches the decline after peak in the optical bands until ∼ 80 days after explosion (the latest photometric point of SN 2020kyg available).However, the r114_d6.0_Z_6000K_0.1Nimodel is too bright at peak by up to half a magnitude in the optical bands and a magnitude or more in the NIR bands, while also declining too slowly after peak in the NIR bands.Although we expected the r114_d6.0_Z_6000K_0.1Nimodel would be too bright compared to SN 2020kyg, the more significant overproduction of flux at NIR wavelengths suggests our adopted remnant temperature may be too low.
To test this we carried out another simulation also based on Model r114_d6.0_Zbut with an increased fixed remnant black body tem-perature of 8000 K and a further reduced 56 Ni mass a twentieth of that predicted by L22.This model did produce light curves which were slightly fainter at peak, up to ∼ 0.2 magnitudes in optical bands and as much as ∼ 0.4 in NIR bands, leading to closer agreement with SN 2020kyg.This model also produced NIR band light curves that matched the decline SN 2020kyg more successfully.However, this came at the expense of significantly reducing the agreement of the model decline in the r, i and z-bands and in general providing a worse match to the band light curves of SN 2020kyg compared to the r114_d6.0_Z_6000K_0.1Nimodel.We also note that as we move to faint deflagration models in which we include the remnant contribution, the models become more sensitive to the choice of remnant temperature.This may be a result of the ejecta becoming optically thin more rapidly for faint models meaning the remnant radiation has a larger contribution to the synthetic observables at earlier times for faint models.
Overall, we have been able to achieve a significant improvement in agreement with the band light curves of SN 2020kyg, particularly in decline, by including the remnant radiation in our simulations.However, there is no obvious way using our approach to produce a simulation based on the standard L22 pure deflagration models with remnant contribution included that can produce good agreement with the light curves of SN 2020kyg across all bands (see Section 4 for further discussion).et al. 2021).Also plotted are the r48_d5.0_Zmodel spectra scaled so their peak flux matches that of Model r48_d5.0_Z_6000K_0.33Ni.Times are relative to B-band peak.The spectra of SN 2019muj have been corrected for red shift and reddening using the values from Barna et al. (2021).The observed spectra shown here were taken from WISeREP (Yaron & Gal-Yam 2012).Srivastav et al. 2022).Also plotted are the r114_d6.0_Zspectra scaled so their peak flux matches that of Model r114_d6.0_Z_6000K_0.1Niwhere required, to allow clearer comparison of spectral features.Times are relative to g-band peak.The spectra of SN 2020kyg have been corrected for red shift and reddening using the values from Srivastav et al. (2022).

Spectra
For the last 3 epochs shown (16.4,24.4 and 53.8 days after g peak) the r114_d6.0_Z_6000K_0.1Nimodel provides a significantly better match to the absolute flux of SN 2020kyg compared to the r114_d6.0_Zmodel which becomes much too faint.Additionally, while the spectra of the r114_d6.0_Z_6000K_0.1Nimodel are also slightly too blue they provide a significantly better match to the overall colours: the r114_d6.0_Zmodel is significantly too blue.Therefore, as was the case for intermediate-luminosity deflagration models, including the remnant radiation for faint deflagration models improves agreement with observed SNe Iax because of a wavelength dependent increase in flux for the r114_d6.0_Z_6000K_0.1Nimodel, producing a colour evolution more similar to SN 2020kyg.
The r114_d6.0_Z_6000K_0.1Nimodel is able to reproduce the strength and location of many of the spectral features of SN 2020kyg for the earlier epochs shown.However, for the latest two epochs shown, while the model is still reproducing some of the spectral features, they are often too weak.This is particularly clear for the last epoch shown at wavelengths redder than ∼ 5500 Å where the r114_d6.0_Z_6000K_0.1Nimodel spectra struggles to match the large number of strong sharp features exhibited by the spectrum of SN 2020kyg.This is likely because the spectra of SN 2020kyg at these later epochs show an increasing number of features due to forbidden transitions (Srivastav et al. 2022) which can not be treated accurately using the atomic data set and approximate NLTE treatment we have adopted in our radiative transfer simulations for this work.Adopting the full NLTE treatment and more extensive atomic dataset of Shingles et al. (2020) may produce model spectra which more successfully reproduce these spectral features at later times.
Focusing again on the latest two epochs shown we see interesting differences in the evolution of the Ca ii NIR triplet.In the spectral comparison at 24.4 days after peak the absorption feature at ∼ 8400 Å, which we have confirmed to be the most blue-ward absorption feature of the Ca ii NIR triplet in the models, matches the velocity of the feature in SN 2020kyg well.However, for the spectral comparison at 53.8 days after peak while this absorption feature has remained at the same velocity for the model it has shifted significantly to the red for SN 2020kyg.Additionally, comparing the overall profile of the Ca ii NIR triplet in the latest epoch shown we see that it is both bluer and broader for the deflagration models compared to SN 2020kyg.This explains the lack of emission from the 8542 Å component of the Ca ii NIR triplet for the models as it is suppressed by the blue shifted absorption tail of the 8662 Å component, something which Maeda & Kawabata (2022) also observed for their spectral models of SN 2019muj.These comparisons of the Ca ii NIR triplet suggest there is a low velocity, higher density region of the ejecta in SN 2020kyg, where spectral features are forming at later epochs, which is not accounted for in our models.Such an ejecta structure has already been proposed for other observed SNe Iax (see e.g.Sahu et al. 2008) while Foley et al. (2016) and Shen & Schwab (2017) suggest such an ejecta structure may be consistent with radioactively driven post SNe Iax remnant winds.Alternatively, Maeda & Kawabata (2022) have suggested that a scenario involving a secondary higher density ejecta component which is ejected from the remnant envelope ∼ a month after explosion may explain this spectral behaviour at later times.

Summary
We have presented a new set of radiative transfer simulations for models based on the  Ch CO WD pure deflagration models of L22 in which we also include energy injection from a luminous remnant, as predicted in the explosion models and possibly observed for the SNe Iax, SN 2008ha (Foley et al. 2014), SN 2012Z (McCully et al. 2022), SN 2014dt (Kawabata et al. 2018) and SN 2019muj (Kawabata et al. 2021).While our models require an assumption about the SED of photons from the remnant we found no significant differences in the synthetic observables of the models where we assumed a black body SED with either constant radius or temperature.Therefore, our models are not overly sensitive to the specific choice of effective remnant temperature at a given time.Moreover, our models which produced best agreement with observed SNe Iax all adopted remnant temperatures consistent with late time temperature estimates of observed SNe Iax.We note, however, that if we were to extend our simulations to later times, the remnant temperature predicted for a model with fixed remnant radius will increasingly diverge from a fixed remnant temperature making this choice more relevant.Additionally, at later times the ejecta become significantly more optically thin meaning radiation emitted from the remnant will be less impacted by interaction with ejecta material and thus will likely become more sensitive to the remnant temperature treatment utilised.
Including the remnant 56 Ni mass predicted for the brightest pure deflagration model from the L22 sequence increases the peak bolometric magnitude of the model by ∼ half a magnitude.This means the brightest viewing angles of this model can match the peak bolometric luminosity of the brightest SNe Iax within their uncertainties, something not previously possible for the L22 single-spark pure deflagration models.Additionally, the light curves are in good agreement with the bright SNe Iax, SN 2005hk in peak magnitude, rise and decline in all bands.This represents a significant improvement as previous bright deflagration models presented by K13, K14 and L22 exhibited declines systematically too rapid compared to bright SNe Iax in the red optical and NIR bands.The inclusion of the remnant also improved spectroscopic agreement with SN 2005hk, primarily due to the better overall flux agreement, although for the latest epoch compared the redder spectrum formed also improved the match with the spectral colours of SN 2005hk.However, including the remnant contribution causes a more pronounced change in ejecta conditions for intermediate-luminosity and faint models leading to bigger improvements in their spectroscopic agreement with observed SNe Iax relative to standard deflagration models.
Of the models investigated in this work, intermediate-luminosity deflagration models with remnant contribution included display the largest improvement in agreement with observed SNe Iax relative to standard deflagration models.Including the remnant produces model light curves in good agreement with the intermediate-luminosity SN Iax, SN 2019muj in terms of peak magnitude, rise and decline in all bands.This again represents a significant improvement compared to standard deflagration models, which show light curves that decline systematically too fast after peak in all bands relative to observed SNe Iax.Additionally, including the remnant contribution leads to significantly improved spectral agreement with SN 2019muj due to the better match with the evolution of the SED.
Including the remnant contribution also significantly improves the agreement with the light curve evolution of the faint SNe Iax, SN 2020kyg.This improvement is primarily due to the slower light curve decline after peak in all bands compared to faint standard deflagration models which show declines that are significantly too fast relative to SN 2020kyg.Additionally, as was the case for intermediateluminosity models, including the remnant contribution leads to significantly better spectral agreement with SN 2020kyg due to the improved match to the evolution of its SED.The inclusion of the remnant does, however, make the model light curves too bright at peak relative to SN 2020kyg in all bands and significantly too bright over the whole simulation time in NIR bands.This is to be expected as the faintest member of the L22 model sequence is already slightly brighter than SN 2020kyg at peak.Therefore, to achieve agreement with the faintest observed SNe Iax for deflagration models with remnant contribution included we require models with an ejecta component ∼ a magnitude or more fainter than the faintest model from the L22 sequence.A future investigation of whether it is possible to produce even fainter pure deflagration models, something L22 suggested can be achieved for single-spark deflagration models, would therefore be of interest.

Implications for the remnant
While our models which produce good agreement with the bright SNe Iax, SN 2005hk, adopt the full remnant 56 Ni mass predicted by L22 our models in best agreement with the intermediate-luminosity SNe Iax, SN 2019muj and faint SNe Iax, SN 2020kyg, had remnant 56 Ni masses one third and one tenth of those predicted by L22 respectively (see Table 2).From Table 1 we see that moving to fainter deflagration models the remnant 56 Ni mass predicted increases relative to the ejecta 56 Ni mass.Adopting the entire remnant 56 Ni mass predicted by L22 and assuming all energy emitted from the 56 Ni decay chain in the remnant contributes to the optical display of the faintest L22 deflagration model (r114_d6.0_Z)produces an event of comparable brightness to SN 2019muj (see Section 3.2.1).L22 suggest it should be possible to produce models fainter than r114_d6.0_Z.However, utilising the entire predicted remnant 56 Ni mass in our intermediate-luminosity and faint deflagration models results in a very large change in the overall properties of the models leading to poor agreement with observed SNe Iax.This suggests not all the energy injected in the remnant from radioactive decays can be used to power the light curves of intermediate-luminosity and faint models.This finding is consistent with what Kromer et al. (2015) found for their faint hybrid CONe WD pure deflagration model they compared to the faint SNe Iax, SN 2008ha.
Models of SNe Iax remnants by Shen & Schwab (2017) that predict some of the radioactive energy released powers a wind in the outer envelope of the remnant provide a potential explanation for this.If a large enough fraction of the energy released in the remnant is used to power such a wind, the energy contributing to the light curve will be significantly reduced.Alternatively, Maeda & Kawabata (2022) have proposed a mechanism in which energy from radioactive decays causes hot 56 Ni rich regions of the remnant to be ejected on the timescale of approximately a month after explosion, again a fraction of the energy injected in the remnant is required to power this process.As we do not take into account either of these mechanisms in our simulations these scenarios are roughly equivalent to having a reduced remnant 56 Ni mass for our models (although time variations in the fraction of deposited energy driving these mechanisms will not be captured).
To estimate the potential impact of remnant winds, we calculated the approximate kinetic energy required to drive winds in the envelope of the remnant for the r48_d5.0_Zand r114_d6.0_Zmodels from the L22 sequence, using the SNe Iax remnant wind models computed by Shen & Schwab (2017) as a guide.The Shen & Schwab (2017) models with greatest remnant masses have core and envelope masses of 0.9 and 0.1 M ⊙ respectively, giving slightly lower but similar total remnant masses to the L22 remnant masses predicted for the models used in this study (1.16 to 1.38 M ⊙ ).Although we are unable to differentiate the remnant core and envelope from the L22 simulations we do have approximate yields of the IGEs and IMEs in the remnant which should primarily be found in the envelope and make up a sizeable fraction of the total envelope mass.The IGEs and IMEs predicted for the remnants sum to 0.074 and 0.053 M ⊙ for the r48_d5.0_Zand r114_d6.0_Zmodels respectively.As some of the envelope will be unburnt CO material, which we can not determine the mass of, an envelope mass of 0.1 M ⊙ seems reasonable for these models and we utilise this envelope mass as an initial guess for the mass driven by the radioactive wind in this calculation.The Shen & Schwab (2017) models have average winds speeds up to ∼ 1100 km s −1 which we adopt for our calculation here to obtain an estimated upper limit on the energy used to drive winds in the remnant envelope.
The total energy deposited by the 56  From the calculation discussed above we estimate the total kinetic energy required to drive the remnant wind over the simulation time is 1.20 × 10 48 erg.For both models, this accounts for just over a third of the discrepancy between the energy we inject in the model remnants and the energy we predict would be deposited in the remnants from the L22 remnant 56 Ni mass estimates.For envelope masses higher than the 0.1 M ⊙ we estimated for these calculations an even greater proportion of the energy injected from radioactive decays in the remnant could be driving winds.It is therefore plausible that energy used to drive such winds can explain why only a fraction of the total energy injected in the remnants from radioactive decays is required to power our r48_d5.0_Z_6000K_0.33Niand r114_d6.0_Z_6000K_0.1Nimodel light curves.
Shen & Schwab (2017) also suggest the hot relatively high density regime of the remnant may be so highly ionised that 56 Ni and 56 Co have radioactive decay rates lower than those we typically use in our simulations, which are appropriate for atomic or low ionisation species.This reduction in decay rate is expected to be relevant in highly ionised regimes with densities ≲ 10 5 g cm −3 (see Shen & Schwab 2017 for details).The remnants predicted by the pure deflagration simulations of L22 have maximum central densities of ∼ 10 5 g cm −3 .Therefore, it is plausible that regions of the remnant may inhabit this regime of reduced decay rates providing another possible explanation for our requirement to adopt lower remnant 56 Ni masses than those predicted by L22.Maeda & Kawabata (2022) alternatively suggest that diffusion timescales in the remnant are long enough that they lead to significant delays in the radiation contributing to the optical display of SNe Iax.In their proposed scenario 56 Ni rich material is ejected from the remnant around a month after explosion.They estimate diffusion in this material will become efficient by around 20 days since explosion allowing radiation originating there to begin to contribute to the optical display.However, they suggest it could be ∼ a year before radiation produced from the decay of 56 Ni and 56 Co in the remaining bound material can diffuse efficiently and thus contribute to the optical display of the SNe Iax.As we assume no significant diffusion time in the remnant in our models, this scenario could also explain the reduced remnant 56 Ni masses we must adopt for intermediate luminosity and faint models.
As discussed above, adopting reduced 56 Ni masses in our simulations can provide a rough approximation of the behaviour expected in our models at early times if delayed radioactive decays (Shen & Schwab 2017) or significant diffusion times in the remnant (Maeda & Kawabata 2022) are considered.However, relative to these cases, reducing the remnant 56 Ni mass will underestimate the total radiation emitted by the remnant over time.Additionally, the contribution of the remnant to the optical display at later times will be significantly underestimated.Determining how these processes would influence the synthetic observables predicted by our simulations requires detailed modelling of the remnant structure and evolution, which is outside the scope of this study.However, we speculate that if we were to account for delayed radioactive decays or significant diffusion times in the remnant in our simulations, the remnant would likely have a reduced impact on the model light curves around peak but a greater influence on the light curve decline at later times relative to the simulations presented here.
We finally note that it is also possible that the L22 simulations overestimate the remnant 56 Ni masses.For fainter deflagration models, the mass of burnt material relative to the total WD mass is much smaller, meaning the yields of the burnt material are more uncertain for these models than for brighter more energetic models which burn more material.However, these uncertainties would have to account for very large reductions in the remnant 56 Ni masses predicted by L22 for the models (up to a factor of 10) and it is therefore highly unlikely these uncertainties are relevant here.

Future work
Including the energy injected from a luminous remnant in our radiative transfer simulations leads to significantly improved agreement with the spectra and light curves of observed SNe Iax, strengthening the case of pure deflagration models for SNe Iax.However, outstanding questions remain.Firstly, it is unclear how the assumptions of the radiative transfer impact our results.In this work, we have adopted the NLTE approximation described by Kromer & Sim (2009).However, including a full NLTE treatment of the plasma conditions has been shown to have important effects on the synthetic observables produced by radiative transfer simulations of SNe Ia, especially at later times (Blondin et al. 2013;Dessart et al. 2014;Shingles et al. 2020Shingles et al. , 2022;;Shen et al. 2021;Collins et al. 2023b).In future work we will utilise the full NLTE treatment of Shingles et al. (2020) to investigate the importance of NLTE effects.In particular, adopting the full NLTE treatment will enable us to extend our simulations to the late phases where the remnant contribution becomes more prominent as the ejecta become increasingly optically thin.This will allow us to investigate if the remnant is required to explain the peculiar late time spectra of observed SNe Iax which never become truly nebular, as suggested by Foley et al. (2016).
Important questions also remain around the nature of the remnant.Hydrodynamic explosion simulations of pure deflagrations following the remnants more closely would provide detailed remnant structures and compositions which would improve our understanding of the remnant structure providing insight into its likely evolution and emission.This is highly relevant to determining if the energy injection we utilise in this study is reasonable, particularly for our intermediate-luminosity and faint models that require remnant 56 Ni masses lower than those predicted by L22.As such, it would also be useful to further investigate the possibility of radioactively driven winds and delayed radioactive decays as suggested by Shen & Schwab (2017) and the second ejecta component proposed to be ejected from the remnant envelope by Maeda & Kawabata (2022), in particular focusing on how introducing such mechanisms impacts the overall properties of pure deflagration models and their agreement with observed SNe Iax.A1 for the Bessel UBVRĲHK bands.Note, the lower relative flux of the models in the JHK-bands meant reliable values could not be determined for the line-of-sight dependent light curve quantities due to the significant Monte Carlo noise.Therefore only the angle averaged values are quoted for these bands.Additionally, we do not include any JHK-band light curve quantities for Model r10_d4.0_Z_R_2000Kdue to its peculiar light curve behaviour in these bands which we believe to be unphysical.

Figure 1 .
Figure 1.Angle averaged bolometric light curves for the r10_d4.0_ZL22 model and our new r10_d4.0_Z_R_8000Ksimulation.We also show the 56 Ni decay energy deposition in the remnant (cyan), and the light curve obtained by adding this to the r10_d4.0_Zmodel light curve (red).The bolometric light curve of the bright SN Iax, SN 2005hk (Phillips et al. 2007) is shown for comparison.

Figure 2 .
Figure2.Angle averaged bolometric and UBVRĲH-band light curves for the r10_d4.0_ZL22 model and four new models from this work in which we also treat energy injected by radioactive decays in a luminous remnant in our radiative transfer simulations.Included for comparison are light curves of SN 2005hk(Phillips et al. 2007), a bright SN Iax.As Model r10_d4.0_Z_R_2000Kdid not provide a well converged temperature solution at later times we only plot the light curves for this model up to 65 days.This poor convergence is a result of the low remnant temperature adopted for this model.

Figure 3 .
Figure 3. Absolute flux spectral comparisons between the angle averaged r10_d4.0_Z_R_8000Kmodel and r10_d4.0_ZL22 model along with SN 2005hk(Phillips et al. 2007).Also plotted are the r10_d4.0_Zspectra scaled so their peak flux matches that of the r10_d4.0_Z_R_8000Kmodel spectra.Times are relative to B-band peak.The spectra of SN 2005hk have been corrected for reddening using the values fromPhillips et al. (2007) and redshift corrected using the values fromSahu et al. (2008).

Figure
Figure 6 shows spectroscopic comparisons between Model r114_d6.0_Z_6000K_0.1Ni,Model r114_d6.0_Z,Model r114_d6.0_Zscaled to match Model r114_d6.0_Z_6000K_0.1Niat peak and SN 2020kyg.For the earliest epoch shown (3.1 days before g peak) both models are spectroscopically very similar and provide a good match in absolute flux to the spectrum of SN 2020kyg.At the next epoch shown (3.5 days after peak) the r114_d6.0_Zmodel still shows a good match to the absolute flux of SN 2020kyg while the r114_d6.0_Z_6000K_0.1Nimodel becomes too bright.

Figure 5 .
Figure5.Absolute flux spectral comparisons between the angle averaged r48_d5.0_Z_6000K_0.33Nimodel, r48_d5.0_ZL22 model and SN 2019muj(Barna et al. 2021).Also plotted are the r48_d5.0_Zmodel spectra scaled so their peak flux matches that of Model r48_d5.0_Z_6000K_0.33Ni.Times are relative to B-band peak.The spectra of SN 2019muj have been corrected for red shift and reddening using the values fromBarna et al. (2021).The observed spectra shown here were taken from WISeREP(Yaron & Gal-Yam 2012).

Figure 6 .
Figure6.Absolute flux spectral comparisons between the angle averaged r114_d6.0_Z_6000K_0.1Nimodel, r114_d6.0_ZL22 model and SN 2020kyg(Srivastav et al. 2022).Also plotted are the r114_d6.0_Zspectra scaled so their peak flux matches that of Model r114_d6.0_Z_6000K_0.1Niwhere required, to allow clearer comparison of spectral features.Times are relative to g-band peak.The spectra of SN 2020kyg have been corrected for red shift and reddening using the values fromSrivastav et al. (2022).

Table 1 .
As we move to fainter deflagration models the remnant 56 Ni mass increases relative to the ejecta 56 Ni mass.Therefore, if all energy released from

Table 1 .
Summary of the total and 56 Ni masses for the ejecta and remnant material of the L22 models which our new models are based on.Masses are given in solar masses.

Table A2 .
Same as Table