© The Authors 2023

1. Introduction

Current models for Type Ia supernovae (SNe Ia) invoke variations in the mass of the exploding carbon-oxygen white dwarf (WD) and in the conditions of the thermonuclear runaway. These models include delayed detonations in near-Chandrasekhar-mass (M_Ch) WDs (Khokhlov 1991), double detonations in sub-M_Ch WDs (e.g. Woosley & Weaver 1994), and violent mergers of two sub-M_Ch WDs (e.g. Pakmor et al. 2012). However, due to numerous degeneracies in SN Ia light curves and spectra, distinguishing between these various models has been a challenge (e.g. Maoz et al. 2014).

The difficulty arises in part from the multi-dimensional nature of the explosion while most radiative-transfer simulations assume spherical symmetry (see e.g. Gamezo et al. 2005; Seitenzahl et al. 2016; Raskin et al. 2009; and Pakmor et al. 2010 for examples of 3D explosion models). Another difficulty resides in the intrinsic complexity of radiative transfer in SN Ia ejecta, including non-local thermodynamic equilibrium (non-LTE) simulations and non-thermal effects, or limitations of the atomic data (see e.g. Höflich et al. 1998; Pinto & Eastman 2000; Sim 2007; Dessart et al. 2014; Shen et al. 2021; and Blondin et al. 2022a). These issues are composition dependent and related to the complicated explosive nucleosynthesis in SN Ia explosions (see e.g. Bravo & Martínez-Pinedo 2012; Bravo 2020; and Seitenzahl & Townsley 2017 for a review).

While the early high-brightness phase of SNe Ia (≲50 d post explosion) yields constraints on the ejecta mass and kinetic energy, as well as the yields of intermediate-mass elements (IMEs) and ⁵⁶Ni, the late nebular phase (> 100 d post explosion) can provide complementary information. At such times, the ejecta is powered by ⁵⁶Co decay with an increasing fraction of the deposited energy arising from positrons, in conditions that are close to steady state. The low ejecta density (electron densities n_e ≲ 10⁶ cm⁻³) is conducive to forbidden-line emission which controls the cooling rate. This cooling is dominated by a few strong lines that provide important nucleosynthetic information.

In particular, the abundance of stable iron-group elements (IGE) that are synthesised during the explosion has been proposed to distinguish sub-M_Ch from M_Ch progenitors (e.g. the Ni/Fe ratio; Flörs et al. 2020). Accurate abundance determinations are however difficult to obtain given the amount of line overlap and the sensitivity of line strengths to small variations in the ionisation state (Blondin et al. 2022b). The advent of JWST has paved the way for systematic studies of SNe Ia in the mid-infrared (MIR), where lines are typically less blended (and less sensitive to the electron temperature) than in the optical and near-infrared (NIR). This allows for a more secure identification of lines as well as detailed studies of their morphology, which may constrain the spatial distribution of the emitting material (e.g. Gerardy et al. 2007).

The Type Ia SN 2021aefx was recently observed with JWST using the Near Infrared Spectrograph (NIRSpec) in the Fixed Slits (FS) spectroscopy mode and the Mid Infrared Instrument (MIRI) in the Low Resolution Spectroscopy (LRS) mode (Kwok et al. 2023a). It was discovered within a few hours post explosion in NGC 1566, located at an estimated distance of 18 ± 2 Mpc (and recession velocity cz = 1504 km s⁻¹), with a Galactic reddening E(B​ − ​V)≈0.008 mag (Schlafly & Finkbeiner 2011) and an estimated host-galaxy reddening E(B​ − ​V) = 0.097 mag (Hosseinzadeh et al. 2022). Early-time observations revealed a UV flux excess within the first few days after explosion in this otherwise normal event (Ashall et al. 2022; Hosseinzadeh et al. 2022). The nebular JWST spectrum at +255 d past maximum light (∼270 d post explosion) published by Kwok et al. (2023a) covered the full 0.35–14 μm range and revealed the presence of numerous forbidden lines of IGEs and a distinct line due to argon ([Ar III] 8.99 μm), whose flat-topped profile was interpreted as indicative of a chemically stratified ejecta. Subsequent modelling by DerKacy et al. (2023) of a MIR spectrum taken ∼2 months later showed satisfactory agreement with a near-M_Ch delayed-detonation model with ∼0.6 M_⊙ of ⁵⁶Ni. However, their study did not confront these data to alternative SN Ia explosion models.

Here we compare a diverse set of state-of-the-art, public SN Ia explosion models with the full optical-MIR nebular spectrum of SN 2021aefx, and evaluate which one, if any, fares better. In the next section, we present our numerical approach. We then present our theoretical spectra for each SN Ia explosion model in Sect. 3, and study the impact of ejecta asymmetries in Sect. 4. We discuss the implications of our results in Sect. 5 and conclude in Sect. 6. All model outputs are publicly available online¹.

2. Numerical approach

We select previously-published models from the Heidelberg Supernova Model Archive (HESMA; Kromer et al. 2017) with a ⁵⁶Ni yield in the range 0.5–0.8 M_⊙ (see Table 1), compatible with the amount inferred for SN 2021aefx². These include 17 models but for only four explosion scenarios, namely: a M_Ch delayed detonation, M_Ch pulsationally assisted gravitationally-confined detonations, sub-M_Ch double detonations, and a violent merger of two sub-M_Ch WDs. In the following section (Sect. 3), we select one representative model for each scenario, namely the M_Ch delayed-detonation model N100 of Seitenzahl et al. (2013; HESMA name ddt_2013_N100, hereafter DDT), the M_Ch pulsationally assisted gravitationally-confined detonation model r10_d1.0 of Lach et al. (2022; HESMA name gcd_2021_r10_d1.0, hereafter GCD), the sub-M_Ch double-detonation model M1002 of Gronow et al. (2021; HESMA name doubledet_2021_M1002_1, hereafter DBLEDET), and the violent merger model of Pakmor et al. (2012; HESMA name merger_2012_11+09, hereafter MERGER).

All models correspond to spherically-averaged versions of a 3D simulation. This averaging causes a systematic overestimate of the total mass by 3–8% and a difference of up to ±10% for several isotopic and elemental yields compared to the original 3D models (see Table A.1). For the N100 model we also obtained non-spherically-averaged radial profiles along six directions (I. Seitenzahl, priv. comm.), whose density was rescaled to match the total mass (1.4 M_⊙) of the original 3D model. Here the yields can vary significantly with respect to the original 3D model, by up to a factor ∼2 in some cases (Table A.2). We subsequently discuss these models in Sect. 4.

Starting with the spherically-averaged density and abundance profiles at t ≈ 100 s post explosion available on HESMA, we generated initial conditions at 270 d post explosion taking into account changes in composition induced by the decay of radioactive isotopes (mainly ⁵⁶Co decay at this time) and the decrease in density due to homologous expansion (ρ ∝ 1/t³). We applied a small radial mixing to the HESMA inputs with a characteristic velocity width Δv_mix = 300 or 400 km s⁻¹ to smooth sharp variations in composition (see e.g. Blondin et al. 2022a). The impact of this mixing on the total yields is below ∼0.1% for all important species. The initial temperature was set to 5000 K throughout the ejecta (this value was found to ease the convergence of the radiative-transfer calculation and is not too far from the final value in the inner ejecta layers; see Sect. 3). The density and abundance profiles at 270 d post explosion for a subset of four models are shown in Figs. 1 and 2.

Fig. 1. Density profiles at 270 d post explosion for our reference model set: the delayed-detonation model ddt_2013_N100 (DDT), the pulsationally assisted gravitationally-confined detonation model gcd_2021_r10_d1.0 (GCD), the double detonation model doubledet_2021_M1002_1 (DBLEDET), and the violent merger model merger_2012_11+09 (MERGER).

Fig. 2. Abundance profiles at 270 d post explosion for our reference model set.

We then solve the 1D non-LTE radiative transfer with CMFGEN (Hillier & Dessart 2012) assuming steady state. Non-local energy deposition from radioactive decay (> 99.9% of which results from ⁵⁶Co decay) was treated using a Monte-Carlo approach for γ-ray transport. Non-thermal processes are accounted for through a solution of the Spencer-Fano equation (see Li et al. 2012). A new temperature solution is produced by CMFGEN as part of the full non-LTE solution.

The following ions were included: He I–II, C I–III, N I–III, O I–III, Ne I–III, Na I, Mg II–III, Al II–III, Si II–IV, S II–IV, Ar I–III, Ca II–IV, Sc II–III, Ti II–III, Cr II–IV, Mn II–III, Fe I–V, Co I–IV, and Ni I–V. For all the aforementioned ions, we also considered ionisations to and recombinations from the ground state of the next ionisation stage (e.g. Fe VI in the case of Fe). The number of levels for all ions and a description of the sources of atomic data can be found in Appendix B.

3. Ejecta properties and predicted spectra

3.1. Energy deposition, temperature, and ionisation structure

At nebular times, and specifically at 270 d post explosion in our model set, only ∼5% of the total decay power is deposited in the ejecta, of which about 50–70% comes from positrons, which we assume are absorbed locally (see Penney & Hoeflich 2014 for a discussion on the validity of this assumption). Ejecta layers below ∼16 000 km s⁻¹ capture 99% of the total volume-integrated decay energy deposition rate (with the notable exception of the double-detonation model doubledet_2021_M1010_1, which also has the highest ⁵⁶Ni yield of our model set; see Table 1). This implies that nebular-phase spectra best probe the regions rich in IGEs but provide weaker constraints on the outer ejecta where the bulk of IMEs is located.

All models in our set converge to a roughly uniform temperature of about 6000 K below 10 000 km s⁻¹ (Fig. 3, upper panel), with some diversity at larger velocities that reflect the variations in the outer density (Fig. 1; a lower outer density yields a higher temperature for an equivalent power absorbed). Overall, the higher-mass MERGER model is the cooler one, while the lower-mass DBLEDET model is the hotter one. The ionisation state of IGEs, which dominate the composition, reflect closely this temperature profile (Fig. 3, lower panel). In the MERGER model, Fe is a mixture of Fe⁺ and Fe²⁺, while Fe²⁺ dominates in all other models (see also Fig. C.5). These properties suggest that lines of Fe I or Fe IV should be subdominant. The mean ionisation profiles for Co and Ni (not shown) are qualitatively similar to those of Fe, in part because of their similar ionisation potential. However, Co is typically more ionised, with a near-equal mixture of Co²⁺ and Co³⁺ in the non-MERGER models (Fig. C.6), and the Ni ionisation level is more uniform across all models, with Ni²⁺ being the dominant stage (Fig. C.7).

Fig. 3. Characteristic ejecta properties at 270 d for our reference model set. Top: temperature profiles. Middle: electron density profiles. Bottom: mean Fe ionisation profiles. The mean ionisation is defined as ∑iini+/∑ini+, where ni+ is the number density of ionisation stage i for Fe, such that a mean ionisation of ∼2 indicates that Fe2+ is the dominant stage.

3.2. Model spectra: Overall comparison to SN 2021aefx

Figure 4 shows the synthetic spectra for models DDT, GCD, DBLEDET, and MERGER compared to SN 2021aefx at 270 d post explosion in the wavelength range 0.35–14 μm³. These data have been corrected for redshift and reddening (Kwok et al. 2023a), while the models have been scaled to a distance of 18 Mpc. We report line identifications in Table 2, based on their Sobolev equivalent width (EW)⁴.

Fig. 4. Spectra of our reference model set at 270 days post explosion compared to SN 2021aefx over the wavelength range 0.35–14 μm (logarithmic scale), with flux in Fν on a non-linear (arcsinh) scale, as in Fig. 1 of Kwok et al. (2023a). The SN 2021aefx spectrum has been corrected for redshift and extinction (assuming a host-galaxy reddening of 0.097 mag and a MW reddening of 0.008 mag). The synthetic fluxes correspond to the same assumed distance to SN 2021aefx of 18 Mpc; they have not been rescaled or normalised in any way. We include selected line identifications based on their maximum Sobolev EW (see Table 2). The insets show normalised line profiles in velocity space for selected transitions, illustrating variations in line widths and morphology.

The main lesson from this figure is that the models predict all the observed spectral features of SN 2021aefx from the blue end of the optical until the red end of the JWST range. Unfortunately, no single model produces a perfect match. The modest difference between the spectral properties of these four distinct explosion models reflects the similarity in composition for the dominant coolants in the ejecta. There are, however, numerous offsets in specific lines or features that result from differences in the abundance and ionisation of specific elements. The spherical averaging of the input explosion models will of course impact the predicted spectra. We attempt to qualitatively evaluate the impact of ejecta asymmetries later in Sect. 4.

Given the similar ⁵⁶Ni mass in these models, the deposited decay power and hence the total radiated energy is also similar, and the line emission is strongly coupled. Altering the flux in one line will affect other lines whose flux will adjust to maintain a constant cooling rate. This effect is most important for the strong IGE lines in the optical (namely the [Fe III] complex around 0.5 μm) since they contain most of the flux. The effect is also seen in some IME lines that can be important coolants in the outer IME-rich regions of the ejecta (e.g. [S III–IV], [Ar II–III], and [Ca II, IV]). Conversely, weak lines that contribute a negligible amount to the cooling are not affected by this flux redistribution.

We attempted to quantify the quality of each match by considering various statistical estimators, such as the mean absolute deviation (MAD) or the mean fractional error (MFE) with respect to the observed fluxes over the entire wavelength range (see Appendix D and Table D.1). To minimise the impact of a mismatch in a few lines, we also considered a logarithmic flux scale when computing the MAD or MFE (or the ‘Score’ of Omand & Jerkstrand 2023). Last, we tried switching to flux per unit wavelength as opposed to per unit frequency, since the latter visually enhances the flux in the NIR (1–5 μm) and MIR (5–14 μm) ranges compared to the optical range (0.35–1 μm), even though their combined contribution to the total flux is only ∼6% and ∼7%, respectively (see Fig. 5). However, the results were inconclusive, as each model was identified in turn as the ‘best-fit’ model depending on the adopted estimator or flux scale (also when taking into account the uncertainty in host-galaxy extinction, or after scaling the model fluxes to match the mean observed flux). These statistics therefore only provide a subjective measure of the goodness of fit, and reflect the difficulty in distinguishing a best-match model for SN 2021aefx.

Fig. 5. Normalised cumulative integral of the flux per unit wavelength over the range 0.35–14 μm in SN 2021aefx and in our reference model set. The vertical lines mark the boundaries between the optical-NIR and NIR-MIR ranges as defined in Sect. 3.3.

Throughout the optical and IR ranges, the spectrum is dominated by forbidden lines. Some exhibit clearly a broad flat top, e.g. [Ar III] 8.99 μm (except in the MERGER model) or [Ni II] 6.63 μm (most easily seen in the MERGER model), indicating a formation starting at large velocities and possibly extending to even larger velocities (the corresponding models have little Ar below about 5000 km s⁻¹ or exhibit an ionisation stratification). Numerous lines show a Gaussian-like profile because they form throughout the ejecta (e.g. [Fe III] 0.47 μm). Only a few permitted transitions are predicted in the optical, such as the Ca II NIR triplet (though significantly weaker in the DBLEDET model due to the higher ionisation; see Fig. C.4), Fe II 0.52 μm, or Na I D (only noticeable in the MERGER model).

The ‘flat’-top profiles in the models exhibit a slight blue excess which we have identified to arise from the relativistic boost to the emissivity. When neglecting such a relativistic effect, the line exhibits a slight flux excess in the red, which is then caused by the frequency redistribution associated with electron scattering (see Appendix E). The observed slant in [Ar III] 8.99 μm is biased to the red and is probably dominated by something lacking in our models (e.g. incomplete atomic data or ejecta asymmetry; DerKacy et al. 2023). We note, however, that line overlap can also reverse the direction of this tilt (see Sect. 3.3.3 below).

All lines from IGEs have typical line widths in the range ∼8000–14 000 km s⁻¹ (FWHM values from Kwok et al. 2023a), while those of IMEs, which tend to be present at large velocities only, exhibit significantly larger widths (FWHM ≳ 20 000 km s⁻¹). Nearly all lines appear centred around their rest wavelength, and the only features exhibiting a significant offset ≳2000 km s⁻¹ are affected by line overlap (e.g. [Fe III] 3.23 μm, [Ni II] 6.63 μm, [Ar II] 6.98 μm; see Kwok et al. 2023a), with the notable exception of [Ni III] 7.35 μm (see Sect. 3.3.3 below). Ejecta asymmetry is invoked by DerKacy et al. (2023) but the profile skewness observed, for example, in individual lines such as [Co III] 11.89 μm is much weaker than that identified in multiple lines in nebular phase spectra of some SNe Ia (see e.g. Dong et al. 2015). The inset in Fig. 4 gives the typical profile morphologies in our model spectra, using [Ni IV] 8.40 μm, [Ar III] 8.99 μm, and [Co III] 11.89 μm.

3.3. Detailed comparison in restricted wavelength ranges

Figure 6 shows the spectra in restricted wavelength ranges (optical: 0.35–1 μm; NIR: 1–5 μm; lower MIR: 5–14 μm), including the 14–28 μm range (upper MIR) for which there is no data at that epoch. In what follows we discuss each of these ranges in turn.

Fig. 6. Same as Fig. 4 but for restricted wavelength ranges (from top to bottom): optical (0.35–1 μm), NIR (1–5 μm), lower MIR (5–14 μm) and upper MIR (14–28 μm). The ℱX−Y μm label gives the fraction of the total optical to MIR flux (0.35–14 μm) for SN 2021aefx emitted within the wavelength range of each plot. We include line identifications based on their Sobolev EW (see Table 2). Transitions connected to the ground state are marked with a ‘†’ symbol. For consecutive lines within a multiplet (connected by a horizontal line), we give the gf-weighted mean wavelength of the transitions. Lines that only appear in one model class are labelled with the corresponding colour.

3.3.1. Optical range (0.35–1 μm)

The optical spectrum is dominated by forbidden lines of doubly-ionised elements (i.e. [Fe III], [Co III], [Ar III], [Ni III], and [S III] lines), with line contributions from singly-ionised Fe, Co, and Ca (both forbidden and permitted transitions) and neutral Fe and Na (semi-forbidden and permitted transitions) mostly in the MERGER model due to the lower ionisation. The (permitted) Ca II NIR triplet is predicted in all models, except the DBLEDET model due to its higher ionisation. This model matches best the [S III] feature at 0.9–1 μm, as both the abundance and ionisation stratification conspire to confine the line-formation region to a narrower velocity range (and hence a more restricted volume) compared to the other models.

As with nebular-phase modelling of SN Ia spectra by other groups (e.g. Mazzali et al. 2015; Shingles et al. 2020), none of our models are able to match the flux level in the ∼0.43 μm feature. In all our models, this feature is dominated by lines of singly-ionised Fe (forbidden lines in the non-MERGER models, with contributions from permitted lines in the MERGER model). The top panel of Fig. 6 seems to suggest that the semi-forbidden line of Fe I] 0.44 μm dominates this feature in the MERGER model, however the combined emission of weaker lines of Fe II is more than a factor of three larger than the combined emission from neutral Fe. These weaker lines are not reported in Table 2 because their Sobolev EW is below our adopted EW cut.

3.3.2. NIR range (1–5 μm)

The NIR spectral range (second panel of Fig. 6) is dominated by forbidden transitions of Fe II and Co II and all four models predict the observed features with the notable exception of the weak featureless continuum beyond 4 μm (which could possibly arise from molecule emission or dust formation in the outer ejecta⁵; see e.g. Jerkstrand et al. 2012 in the core-collapse SN context). The MERGER model indicates that a lower ionisation would improve the agreement for the Fe II and Co II lines in the three other models. Higher ionisation features due to [Fe III] 2.24, 3.01 μm, [Ca IV] 3.21 μm and [Ni III] 3.39, 3.80 μm are also present in all models. We note that the [Co III] 1.55 μm line identified by Kwok et al. (2023a) is due to Co II in our models. Moreover, none of our models predict a contribution from [Ni I] 3.12 μm to the 3.2 μm feature (see Sect. 5.3), which is dominated by [Ca IV] 3.21 μm in our models (Nahar & Shafique 2023), with a modest contribution from [Fe III] 3.23 μm.

3.3.3. Lower MIR range (5–14 μm)

The lower MIR range (third panel of Fig. 6) reveals the presence of lines from three distinct ionisation stages (II–IV) as in the NIR range, and lines from the IMEs argon and sulfur as in the optical range. We also tentatively identify the 12.8 μm feature as [Ne II] 12.81 μm, which is particularly strong and centrally peaked in the MERGER model due to the presence of Ne in the inner ejecta layers (see Fig. 2), combined with the low ionisation. The strongest lines in this range are low-lying forbidden transitions either directly connected to the ground state (marked with a ‘†’ in Fig. 6 and Table 2), such as the prominent [Ni III] 7.35 μm and [Co III] 11.89 μm lines, or part of a multiplet including the ground state, such as [Ni III] 11.00 μm. The [Ni IV] and [S IV] transitions connect levels at most ≲0.3 eV (equivalent electron temperature ≲3500 K) above the ground state which are easily collisionally excited given the ejecta temperature in the corresponding layers.

In the non-MERGER models, the sharp drop in the IME abundance below ∼3000–4000 km s⁻¹ (∼8000 km s⁻¹ in the DBLEDET model) results in a broad flat-top morphology for lines of Ca IV, Ar II–III and S IV, while the presence of IGEs down to the centre results in centrally-peaked lines of Fe III, Co II–III, and Ni III–IV. The weak [Ni II] 6.63 μm line is flat-topped in these models owing to an ionisation (as opposed to abundance) stratification. In the MERGER model, the situation is reversed as IMEs dominate the innermost regions while IGEs are present above ∼1000–2000 km s⁻¹, resulting in Gaussian-like profiles for the Ca, Ar and S lines (e.g. [Ar II] 6.98 μm), and narrow flat-top profiles for the Ni, Co, and Fe lines (e.g. [Co III] 11.89 μm). We note that line overlap can reverse the direction of the slant of the flat top (from blue to red) in some profiles, as seen for instance in [Co II] 10.52 μm (overlap with [Ni II] 10.68 μm) or [Ni III] 11.00 μm (overlap with [Co II] 11.16 μm).

The strength of the [Ni III] 7.35 μm line is largely overestimated in the DDT and MERGER models (factor of ∼4 at line centre), as is the [Ni III] 11.00 μm line (factor of ∼2 at line centre). These lines are particularly sensitive to the collisional strengths for Ni III (see Sect. 5.1), which set the collisional de-excitation rate from the upper level and hence the critical density for each transition (n_crit ≈ 1.5 × 10⁶ cm⁻³ and ∼6.6 × 10⁵ cm⁻³ at 5000 K for the 7.35 μm and 11.00 μm transitions, respectively). The DDT model has the largest stable Ni abundance of our reference model set (∼0.08 M_⊙; see Table 1), with a peak in its abundance profile around 9500 km s⁻¹ (Fig. 2, top panel) where the electron density drops below n_crit (Fig. 3, middle panel). The stable Ni abundance is lower in the MERGER model (∼0.03 M_⊙) and the electron density is higher, but the offset of the IGE distribution to larger velocities ensures there is a sufficient amount of (doubly-ionised) Ni in layers where n_e ≲ n_crit, resulting in strong [Ni III] lines.

The [Ni III] 7.35 μm line also displays one of the largest kinematic offsets (3000 ± 1400 km s⁻¹; Kwok et al. 2023a) that cannot be attributed to line overlap, and which our 1D modelling approach cannot accommodate. However, Kwok et al. (2023a) note that the MIRI wavelength solution is more uncertain in this lower-wavelength range, and the measured offset for the redder isolated [Ni IV] 8.40 μm line is consistent with no offset (1300 ± 1200 km s⁻¹).

3.3.4. Upper MIR range (14–28 μm)

The upper MIR range (fourth panel of Fig. 6) is shown for comparison with future JWST observations extending to higher wavelengths. The key spectroscopic signatures comprise a double-hump feature centred on 16 μm due to [Co IV] 15.64 μm (and [Ne III] 15.55 μm in the MERGER model) and [Co III] 16.39 μm, a broad feature around 18–19 μm dominated by [S III] 18.71 μm with a contribution on the blue side from [Fe II] 17.93 μm (whose blue-to-red tilt in the MERGER model is due to line overlap), and a prominent line at ∼23 μm due to [Fe III] 22.92 μm (with a small contribution due to [Co IV] 22.79 μm). This [Fe III] line constitutes almost 50% of the total flux in this range for the non-MERGER models (and ∼30% for the MERGER model). It was speculated to account for most of the observed flux in the F2100W filter (spanning roughly 18.5–23.5 μm) in JWST imaging of SN 2021aefx taken at +357 d past maximum by Chen et al. (2023). In our model spectra corresponding to ∼100 days earlier, this line constitutes 75–80% of the flux in the F2100W bandpass. Additional lines of Co II and Fe II contribute to the spectrum of the MERGER model due to the lower ionisation. The analysis of the differences in line-profile morphology in the lower MIR range (5–14 μm; Sect. 3.3.3) also applies here.

4. Predicted spectra for different rays of a 3D model

One fundamental issue with spherical averaging of multi-dimensional explosion models is the microscopic mixing that results in otherwise chemically-segregated zones at a given radius but different ejecta-centred directions (see e.g. Dessart & Hillier 2020 in the context of core-collapse SNe). In SNe Ia with significant ejecta asymmetries, with WD-WD mergers and collisions being extreme cases, such mixing can alter the composition in the ⁵⁶Co-rich zones at nebular epochs, where local energy deposition by positrons largely governs the plasma emissivity. We attempted to qualitatively evaluate the impact of such asymmetries on our 1D modelling approach by considering different directions in the original 3D version of the DDT model (model N100 in Seitenzahl et al. 2013; I. Seitenzahl, priv. comm.). Specifically we considered the three orthogonal axes of their Cartesian grid $(\widehat{\mathbf{x}},\widehat{\mathbf{y}},\widehat{\mathbf{z}})$ in both positive and negative directions, resulting in six distinct radial profiles.

The density profiles resulted in total ejecta masses in the range 1.15–1.52 M_⊙ (see Table A.2), so we rescaled the density in each case to match the total mass of the original 3D model (1.40 M_⊙). The resulting ⁵⁶Ni mass is in the range 0.50–0.70 M_⊙, corresponding to a ±0.1 M_⊙ (±17%) difference compared to the original 3D model (in which M(⁵⁶Ni) ≈ 0.60 M_⊙; Seitenzahl et al. 2013). As seen in Table A.2, The impact on the yields of specific isotopes and elements is occasionally larger than 50%, as is the case for ⁵⁸Ni in the model corresponding to the positive $\widehat{\mathbf{z}}$ direction (+81%), or Ca in the model corresponding to the positive $\widehat{\mathbf{y}}$ direction (+87%).

We show the density profiles for these six directions in Fig. 7. The largest deviations from the spherically-averaged model occur in the positive y-direction, with a ∼70% higher density at ∼1000 km s⁻¹ and a ∼35% lower density at ∼9000 km s⁻¹. As expected, the deviation in composition from the spherically-averaged model at a given velocity can be far more significant (Fig. 8). While the mass fraction of IGEs (Ni, Co, and Fe) are within ∼50% of one another in the inner ejecta (v < 4000 km s⁻¹), variations can span a factor of a few out to ∼12 000 km s⁻¹ and orders of magnitude farther out. For the IMEs (illustrated here using Ca, Ar, and S) the relative variation in mass fraction at a given velocity is systematically larger, in part due to their underabundance compared to IGEs.

Fig. 7. Density profiles for the spherically-averaged ddt_2013_N100 (DDT) model and along the three orthogonal axes of the original 3D Cartesian grid, in both positive ({x,y,z}pos) and negative ({x,y,z}neg) directions. The bottom panel shows the density ratio with respect to the DDT model.

The variation in spectral properties along the six directions compared to the spherically-averaged model partly reflects these variations in abundance (Fig. 9). For instance the width of the prominent [Ar III] 8.99 μm line is set by the extent of the Ar hole (which determines the width of the flat-top component⁶; see Sect. 5.5) and the peak or global offset of the Ar distribution (which determines its FWHM). Thus the model corresponding to the positive $\widehat{\mathbf{y}}$ direction has both the broadest flat-top (approximately ±3000 km s⁻¹ about line centre) and the largest FWHM ≈ 21 200 km s⁻¹ (see inset (2) in Fig. 9). Conversely, the model corresponding to the positive $\widehat{\mathbf{z}}$ direction displays a flat-top extending only ±2000 km s⁻¹ about line centre and a FWHM ≈ 13 600 km s⁻¹. Modulations in the Ar abundance distribution also leave an imprint on the line profiles: the dip around ∼8000 km s⁻¹ in the positive $\widehat{\mathbf{x}}$ direction causes a ‘shoulder’ in the [Ar III] 8.99 μm line profile around ±8000 km s⁻¹ in this model.

The remaining variation results from differences in the temperature and ionisation structures due to variations in the ⁵⁶Co distribution along the different directions, in particular above ∼4000 km s⁻¹ (see upper-right panel in Fig. 8)⁷. The impact is most readily seen in the optical Fe III-dominated complex around ∼0.5 μm. It is strongest in the models corresponding to the positive $\widehat{\mathbf{x}}$ and $\widehat{\mathbf{y}}$ directions, where the Fe ionisation is the largest in the range ∼5000–10 000 km s⁻¹ due to the higher ⁵⁶Co abundance in those layers. This ionisation effect is further enhanced by the larger abundance of Fe in these same layers, which follows that of Co. By comparing Figs. 4 and 9 we see that the spectroscopic variation along different lines of sight in the DDT model can be as large as when considering different explosion models. This impacts the accuracy of abundance determinations when using spherically-averaged ejecta as well as our ability to make quantitative comparisons with observations.

Fig. 8. Abundance profiles for Ni, Co, Fe, Ca, Ar, and S in the spherically-averaged ddt_2013_N100 (DDT) model and along the three orthogonal axes of the original 3D Cartesian grid, in both positive ({x,y,z}pos) and negative ({x,y,z}neg) directions. The lower panels show the logarithm of the abundance ratio with respect to the DDT model.

Fig. 9. Similar to Fig. 4 but for the ddt_2013_N100 (DDT) model along the three orthogonal axes of the original 3D Cartesian grid, in both positive ({x,y,z}pos) and negative ({x,y,z}neg) directions. We also show the original (spherically-averaged) DDT model.

While our approach captures broad trends associated with ejecta asymmetries, it is inherently 1D and the line-of-sight dependence is likely overestimated. Moreover, all the predicted line profiles remain largely symmetric and present no intrinsic Doppler shifts, as inferred for several lines by Kwok et al. (2023a). However, barring a few outliers strongly affected by line overlap, most lines display a kinematic offset either consistent with zero or ≲1000 km s⁻¹ (Table 2 in Kwok et al. 2023a). Uncertainties in the wavelength calibration and the low resolution of the MIRI/LRS spectrograph dominate the measurement uncertainty. For instance, the strong [Co III] 11.89 μm line is offset by 500 ± 900 km s⁻¹ in the +255 d past maximum spectrum analysed here according to Kwok et al. (2023a), and by 740 ± 544 km s⁻¹ in the +323 d past maximum spectrum obtained by DerKacy et al. (2023). This latter measurement, only ∼1.4σ away from a zero velocity offset, is used as evidence for an off-centre distribution of ⁵⁶Ni in the ejecta⁸. Conversely, their model (as does ours) fails to reproduce the much larger offset observed in [Ni III] 7.35 μm (see Sect. 3.3.3), as the stable Ni is concentrated in the centre of their off-centre delayed-detonation model. As noted by the authors of that study, higher-resolution spectra (with the MIRI medium-resolution spectrometer, or MRS) and a more accurate wavelength calibration are needed to reliably measure line shifts in nebular SN Ia spectra taken with JWST. If confirmed for multiple Ni lines, such a large kinematic offset would indicate a large asymmetry in the stable Ni distribution and provide important clues on the explosion mechanism (see e.g. Maeda et al. 2010).

5. Discussion

In what follows we discuss in more detail the impact of uncertainties in the atomic data, focusing on the Ni III collision strengths which greatly affect the lines at 7.35 μm and 11.00 μm (Sect. 5.1). We present a case for the tentative association of the feature at ∼12.8 μm with [Ne II] 12.81 μm in Sect. 5.2. In Sect. 5.3 we discuss the notable absence of lines from neutral ions in our model spectra and whether a lower ionisation through clumping can produce a better match to the observations. Last, we discuss the constraints on the ejecta (and hence progenitor) mass that can be inferred from the strength of stable Ni lines (Sect. 5.4) or the width of the [Ar III] 8.99 μm line (Sect. 5.5).

5.1. Uncertainties in atomic data: Ni III collision strengths

Accurate atomic data are of paramount importance to making reliable predictions with radiative-transfer simulations. The quality of such data is highly variable, and in some instances no data is available. This is particularly true of bound-bound collisional cross-sections between low-lying levels of some ions, relevant to the formation of forbidden lines in the NIR and MIR ranges. For the present calculations, we updated the collisional data used by CMFGEN for the following ions: C I, O I, Ar II, Ca IV, Ti II–III, Co I, Co III, and Ni III–IV.

Since there are no published collisional strengths for transitions within the lowest LS term of Ni III (the 3d8 ³F configuration, which gives rise to the 7.35 μm and 11.00 μm forbidden lines), we carried out calculations of collisional cross-sections for electron-impact excitation of Ni III following the methods of Ramsbottom et al. (2007) and Storey et al. (2016). We report effective collision strengths at 5000 K for transitions among the three lowest levels in Table 3 (a more complete set of values is shown in Table F.2). We also include values based on the calculations of Bautista (2001; M. Bautista, priv. comm.)⁹. As shown in Fig. 10, the predicted line strengths for these three independent sets of collision strengths agree very well with one another. In Appendix F we compare our values with the approximations of Axelrod (1980) for forbidden lines.

Fig. 10. Impact of different values of the effective collision strengths amongst the lowest 3F term in Ni III (see Table 3) on the lines at 7.35 μm and 11.00 μm in our reference model set. In the top row we show the DDT and MERGER models where these lines are largely overestimated. In the bottom row we show the GCD and DBLEDET models where the predicted strengths for both lines are comparable to those seen in SN 2021aefx (note the smaller ordinate range compared to the top row).

As noted previously in Sect. 3.2, collisional excitation of [Ni III] results in a significant overestimate in the strength of the two IR lines at 7.35 μm and 11.00 μm for the DDT and MERGER models. While this visually degrades the match to the SN 2021aefx observations, this still enables a relative comparison of the different models. To illustrate the impact of collisional excitation on the strengths of both [Ni III] lines, we recomputed the spectra for our reference model set by fixing the collision strengths to a constant Υ_lu = 0.1 within the ground configuration (yellow line in Fig. 10). The flux in the [Ni III] lines decreases significantly, by a factor ∼4–7 for the 7.35 μm line and ∼3–4 for the 11.00 μm line. The match of the DDT and MERGER models to SN 2021aefx is actually improved with these artificially low collision strengths, but is significantly degraded for the GCD and DBLEDET models, for which our computed values yielded [Ni III] line strengths comparable to those seen in SN 2021aefx. One cannot therefore invert the problem and determine collision strengths based on a goodness-of-fit test of the models to the data, as this requires a complete control of the underlying systematic uncertainties in the modelling procedure that we are currently lacking.

5.2. Possible identification of [Ne II] 12.81 μm

The broad feature at ∼12.8 μm is tentatively associated with [Ni II] 12.73 μm by Kwok et al. (2023a). In the MIR spectrum taken ∼2 months later published by DerKacy et al. (2023), the feature has somewhat weakened and is associated with a combination of [Fe II] 12.64 μm, [Co III] 12.68 μm, and the aforementioned [Ni II] 12.73 μm line, though it is difficult to gauge the match of the predicted model flux to the observed spectrum as both are shown normalised to the peak of the [Co III] 11.89 μm line (their Fig. 9). Our DDT and MERGER models also predict a [Ni II] 12.73 μm line, but it is at least one order of magnitude too weak compared to the observed flux. If this feature were entirely due to Ni II, it should be matched by our MERGER model which otherwise reproduces the flux levels of other prominent [Ni II] lines (at 1.94 μm and 6.63 μm). We therefore propose an alternative identification for this feature.

Our MERGER model displays a very strong line due to [Ne II] 12.81 μm, whose peaked profile reflects the presence of neon all the way to the innermost region of the ejecta (Fig. 2). This is a natural expectation of this model as the central ejecta is dominated by the ashes of the secondary WD which burns to O and IMEs (Pakmor et al. 2012). The GCD model also displays a moderately strong [Ne II] line. In this model the Ne mass fraction exceeds 10⁻² above ∼10 000 km s⁻¹ and Ne⁺ is the dominant ionisation stage in the range ∼8000–12 500 km s⁻¹, resulting in a broad (FWHM ≈ 20 200 km s⁻¹) flat-top profile. This is inconsistent with the 12.8 μm feature seen in SN 2021aefx, which is centrally peaked and would suggest that Ne is present in the inner ejecta. By artificially setting the minimum Ne mass fraction to 10⁻² all the way to the centre of the GCD model, we are able to reproduce both the strength and morphology of the 12.8 μm feature in SN 2021aefx (Fig. 11). Both models also predict a line due to [Ne III] 15.50, which dominates the ∼15.5 μm feature in the MERGER model. In the other models the feature results from an overlap of [Co II] 15.46 μm and [Co IV] 15.64 μm.

Fig. 11. Impact of the minimum neon mass fraction on the [Ne II] 12.81 μm line in the GCD model. We can reproduce the observed 12.8 μm feature in SN 2021aefx by artificially setting the minimum Ne mass fraction to 10−2 throughout the ejecta.

The [Ne II] 12.81 μm line was first mentioned in the context of SNe Ia by Gerardy et al. (2007) for the low-luminosity SN 2005df, but it was associated with residual background emission from the host galaxy. There are hints of narrow emission features at 12.7–12.8 μm in the MIR spectrum of SN 2014J at ∼137 d post explosion published by Telesco et al. (2015), but the authors associate them with [Co III] 12.68 μm and [Ni II] 12.73 μm. However, the low signal-to-noise ratio and limited resolution (∼0.15 μm, corresponding to ∼3500 km s⁻¹ at 12.8 μm) prevents a thorough investigation of these features. Interestingly, Telesco et al. (2015) do predict a moderately-strong [Ne II] 12.81 μm line in their reference model (see their Table 2), but their synthetic spectra are featureless at this wavelength (their Fig. 4).

Pure deflagration models are also expected to display a prominent [Ne II] 12.81 μm line, due to the large-scale turbulent mixing during the explosion that results in the presence of significant amounts of O, Ne (and unburnt C) in the innermost layers of the ejecta (see e.g. Fink et al. 2014). However, the lower explosion energy of such models results in significantly narrower line profiles compared to what is seen in SN 2021aefx.

5.3. Absence of lines from neutral ions: Impact of clumping

Apart from the MERGER model which displays lines of Fe I and Na I, none of our other models display lines from neutral ions. In particular, all our models lack the [Ni I] 3.12 μm line invoked by Kwok et al. (2023a) to explain the blue edge of the prominent ∼3.2 μm feature (due to [Ca IV] 3.21 μm with a modest contribution from [Fe III] 3.23 μm in our models), and the [Ni I] 7.51 μm line predicted by DerKacy et al. (2023) to dominate over the neighbouring [Ni III] 7.35 μm line in their model for the +323 d past maximum spectrum¹⁰. DerKacy et al. (2023) also predict weaker forbidden lines due to Co I throughout the NIR-MIR range. For these lines to show up in our models, the ionisation would have to be lowered significantly, below that of the MERGER model in the inner ejecta (≲6000 km s⁻¹, where n(Co⁰⁺)/n(Co) < 10⁻³; Fig. C.6).

A clumpy structure may accommodate a broader range of ionisation than currently predicted in our smooth, quasi-homogeneous ejecta (see e.g. Mazzali et al. 2020; Wilk et al. 2020). We have tested the impact of clumping on the DDT model and find that [Ni I] lines only start to emerge for a volume-filling factor f = 0.1, which results in a ten-fold increase of the density in the clumps (Fig. 12). The strongest [Co I] line at 12.25 μm is hardly noticeable in the red wing of the neighbouring [Co III] 11.89 μm line. The overall fit to the SN 2021aefx spectrum is significantly degraded, as lines from doubly and triply-ionised elements become too weak, and those from singly-ionised elements too strong (although the Fe II-dominated feature at ∼0.43 μm remains too weak for all values of the volume-filling factor; likewise, the Fe II 0.517 μm line starts to dominate over [Fe III] 0.527 μm for f < 0.5, but it contributes at most ∼30% of the flux at that wavelength). This clearly indicates that if clumping is present, it cannot be uniform as assumed here. Moreover, it is unclear which physical mechanism, if any, could lead to such strongly-clumped ejecta in SNe Ia.

Fig. 12. Similar to Fig. 4 but for the ddt_2013_N100 (DDT) model with different values for the volume-filling factor f used to approximate a uniformly-clumped ejecta (see Wilk et al. 2020). We show the original (unclumped) DDT model and clumped models with f = 0.75, 0.5, 0.25, 0.1, and 0.05. We highlight lines of neutral ions that emerge in the most clumped models, in particular a strong [Ni I] 3.12 μm line.

As noted previously in Sect. 3.3, the MERGER model suggests that a lower ionisation would improve the agreement for the Fe II and Co II lines in the three other models. This ‘overionisation problem’ in theoretical nebular SN Ia spectra was investigated in the context of sub-M_Ch models by Shingles et al. (2022). They find that artificially reducing the non-thermal ionisation rates leads to a better agreement with observed spectra, although a physical justification for reducing these rates is yet to be found. We have not investigated this effect in the present study.

5.4. Nickel line strengths and stable Ni abundance

DerKacy et al. (2023) argue that the strength of the lines of Ni seen in SN 2021aefx (necessarily resulting from stable isotopes at this time) are evidence for burning at densities > 5 × 10⁸ g cm⁻³, which in turn implies a WD mass of at least ∼1.2 M_⊙. Their M_Ch model synthesises ∼0.06 M_⊙ of stable Ni and underpredicts the strengths of Ni II–IV lines in the MIR spectrum of SN 2021aefx at +323 d past maximum, which would suggest an even higher stable Ni mass. The M_Ch DDT model studied here synthesises ∼0.08 M_⊙ of stable Ni and also underpredicts the strength of Ni II lines in the +255 d past maximum spectrum of SN 2021aefx, but overpredicts that of Ni III lines (Sect. 5.1), which points to an overestimation of the ionisation (see also Shingles et al. 2022). The MERGER models also displays strong Ni lines despite the factor of ∼2.4 lower stable Ni abundance (∼0.03 M_⊙) compared to the DDT model. This is simply the result of the lower ionisation in the MERGER model, which leads to stronger lines of Ni II (and weaker lines of Ni IV) compared to the DDT model (see also Blondin et al. 2022b).

When considering the full model set, however, there is a strong correlation between the total integrated luminosity in Ni I–IV lines and the stable Ni mass (Fig. 13). The correlation is mostly driven by the [Ni III] lines, which account for ∼80–90% of the total Ni luminosity in our models (see Table 4). The correlation is significantly degraded when considering only lines of [Ni II] (Pearson correlation coefficient r = 0.61; cf. last line in Table 4), and highlights the large uncertainties associated with abundance determinations based on lines from a single ionisation stage.

Fig. 13. Total integrated luminosity in all lines of Ni I–IV in the wavelength range 0.35–14 μm versus stable Ni mass in our complete model set. There is a strong correlation between both quantities, with one notable outlier (the MERGER model). The reference delayed-detonation model (DDT) is highlighted with a grey circle.

The MERGER model remains a notable outlier in Fig. 13, and illustrates the impact of ionisation effects when estimating absolute stable Ni yields. This integrated luminosity is complicated to determine observationally because of line overlap, which affects the [Ni II] 6.63 μm, [Ni III] 7.35 μm, and [Ni III] 11.00 μm lines, among others. Despite its weaker correlation with Ni mass, the [Ni IV] 8.40 μm line is more isolated and could serve as a reliable tracer of the stable Ni abundance, provided the ionisation balance is well constrained.

Relating the stable Ni yield to the mass of the exploding WD is complicated. While M_Ch models tend to synthesise more stable Ni at a given ⁵⁶Ni yield and initial metallicity (Blondin et al. 2022b), the lower central density of the M_Ch pulsationally-assisted gravitationally-confined detonation models of Lach et al. (2022) results in a stable Ni yield comparable to sub-M_Ch models (∼0.02–0.03 M_⊙; see Table 1).

5.5. Width of [Ar III] 8.99 μm line and WD mass

Kwok et al. (2023a) fit a thick shell emission profile to the broad (FWHM ≈ 23 700 km s⁻¹) boxy profile of the [Ar III] 8.99 μm line and constrain the shell to be located between 8700 ± 200 km s⁻¹ and 13 500 ± 300 km s⁻¹. This suggests a quasi-absence of Ar in the inner ∼9000 km s⁻¹ of the ejecta, or ‘Ar hole’, compatible with the near-M_Ch model of DerKacy et al. (2023), who use this as additional evidence in favour of a high-mass progenitor for SN 2021aefx. The authors argue that this Ar hole cannot extend to such large velocities in explosions of lower-mass WDs, based on the sub-M_Ch pure detonation model DET2 of Höflich & Khokhlov (1996) in which this hole only extends to ∼6000 km s⁻¹.

In the DBLEDET model from a ∼1 M_⊙ progenitor which we consider here, however, the Ar is located at larger velocities compared to the M_Ch DDT model (Fig. 14, right panel), resulting in a broader [Ar III] 8.99 μm line (FWHM ≈ 21 800 km s⁻¹ cf. ∼16 700 km s⁻¹ for the DDT model; see left panel in Fig. 14). The same holds for other sub-M_Ch pure detonation and double detonation models for 1 M_⊙ WDs available in the literature (e.g. 1D models of Shen et al. 2018; 2D models of Townsley et al. 2019; 3D models of Gronow et al. 2021). Therefore, sub-M_Ch models cannot be excluded based on the large width of the [Ar III] 8.99 μm line observed in SN 2021aefx.

Fig. 14. Left: [Ar III] 8.99 μm line profiles in the MCh DDT model compared to sub-MCh double-detonation models with varying He-shell masses (the reference DBLEDET model has the lowest He-shell mass of 0.02 M⊙). Right: Ar distribution, highlighting the larger extent of the Ar hole in these sub-MCh models compared to the MCh DDT model.

The representative DBLEDET model for the class of double-detonation models consists of a 1 M_⊙ C-O core and a 0.02 M_⊙ He shell. When considering models with the same core mass but more massive He shells (0.03–0.1 M_⊙), we note that the argon distribution is indeed shifted to higher velocities (Fig. 14). The FWHM of the [Ar III] 8.99 μm line increases from ∼21 800 km s⁻¹ for M_shell = 0.02 M_⊙ to ∼26 700 km s⁻¹ for M_shell = 0.1 M_⊙. The trend is reversed for the intermediate-mass shells of 0.03 M_⊙ and 0.05 M_⊙, but it is worth investigating in a more systematic manner.

6. Conclusions

We have compared four classes of public state-of-the-art SN Ia explosion models to nebular observations of SN 2021aefx covering the full 0.35–14 μm range (Kwok et al. 2023a). The input models include M_Ch delayed detonations, pulsationally assisted gravitationally-confined detonations, sub-M_Ch double detonations and a violent WD-WD merger. They were selected from the public HESMA archive to match the ⁵⁶Ni yields expected for normal SNe Ia (∼0.5–0.8 M_⊙). The spherically-averaged density and abundance profiles served as initial conditions to 1D non-LTE radiative-transfer simulations with CMFGEN at 270 d post explosion.

Our main result is that no single model emerges as an obvious candidate for SN 2021aefx based on these data alone. All models predict the same set of spectroscopic features, all of which have an observed counterpart. Conversely, all models lack specific characteristics of the observed spectrum, such as the overall mismatch in lines of singly-ionised IGEs throughout the optical and infrared. We tentatively associate the feature at ∼12.8 μm with [Ne II] 12.81 μm. If confirmed, this would constitute the first firm identification of this line in SN Ia ejecta, and would suggest that neon is present all the way to the innermost region of the ejecta. We also predict a neon line at longer wavelengths due to [Ne III] 15.55 μm, although it overlaps with stronger neighbouring [Co II] and [Co IV] lines.

Differences in the abundance structure amongst the different models (and for different directions in the 3D DDT model) affects the widths and morphology of some spectral lines, such as the prominent [Ar III] 8.99 μm line. The predicted blue tilt of its flat-top profile is due to a relativistic effect in our models. The observed profile is instead tilted towards the red, indicating a large-scale asymmetry in the argon distribution (DerKacy et al. 2023). We note that line overlap can significantly skew the line profiles (centred at zero Doppler shift in our 1D models), and mimic effects normally attributed to ejecta asymmetries.

The largest variations in our model spectra result from differences in ionisation. A larger density, such as in the inner ejecta of the MERGER model, results in a lower ejecta temperature and ionisation state, further enhanced by the increased recombination rate. Variations in the ⁵⁶Co distribution directly impact the decay energy deposition rate, of which a large fraction is deposited (locally in our models) by positrons. Our MERGER model suggest that a slightly lower ionisation would improve the agreement of the other (non-MERGER) models with the SN 2021aefx spectrum. However, none of our models display lines of neutral ions in the NIR or MIR ranges; they only emerge when a significant amount of clumping is introduced, but the match to the observations is then severely degraded.

We further show that the large width of the [Ar III] 8.99 μm line in SN 2021aefx does not invalidate sub-M_Ch models, contrary to claims made by DerKacy et al. (2023). The double-detonation models from a 1 M_⊙ progenitor we consider here display the largest [Ar III] 8.99 μm lines, reflecting the larger extent of the Ar hole in these lower-mass ejecta. Moreover, while the total integrated luminosity in lines of stable Ni strongly correlates with the stable Ni mass, the MERGER model was found to be a clear outlier due its lower ionisation (see also Blondin et al. 2022a). Provided the ionisation balance is well constrained, the isolated [Ni IV] 8.40 μm line could be used to estimate the stable Ni yield in observed SNe Ia. Connecting this yield to the progenitor mass requires an accurate knowledge of the explosion model, as some M_Ch models (e.g. pulsationally-assisted gravitationally-confined detonations) synthesise similar amounts of stable Ni compared to sub-M_Ch models.

As for all radiative-transfer simulations, our results are affected by uncertainties in the atomic data. We could find no published collisional strengths for low-lying forbidden transitions within the lowest ³F term of Ni III, and present them here for the first time. Their values differ significantly from the commonly-used approximations of Axelrod (1980), and the impact on the strength of the [Ni III] lines at 7.35 μm and 11.00 μm is significant. Thus despite the low ejecta densities involved, collisions largely dictate the excitation level and ionisation state of the plasma (including non-thermal collisional ionisation by Compton-scattered electrons; see Shingles et al. 2022). It is clear that more accurate collisional data are needed to correctly interpret these and future infrared observations of SNe Ia.

Throughout this study we have implicitly assumed that the input hydrodynamical models of the explosion provide accurate initial conditions for our radiative-transfer calculations. It is however possible that part of the mismatch with the observations can be attributed to uncertainties in the explosion models themselves, in particular concerning the predicted nucleosynthetic yields and chemical abundance profiles in the ejecta (see e.g. Bravo 2020).

Inferring reliable constraints on the progenitor mass and explosion mechanism of SNe Ia remains a formidable challenge, even with the extended wavelength coverage and exquisite data quality offered by JWST. Combined with early-time and optical observations, nebular infrared spectra of SNe Ia provide an additional validation criterion when comparing a diverse set of SN Ia models to observations.

Acknowledgments

This paper is dedicated to Prof. Tom Marsh (1961-2022), who was S.B.’s lecturer and Bachelor thesis supervisor at the University of Southampton during 1998–2002. S.B. acknowledges support from the Alexander von Humboldt foundation and thanks Sherry Suyu and her group at the Technische Universität München (TUM) for their hospitality. We acknowledge useful discussions with Chris Ashall, Andreas Floers, Peter Hoeflich, Griffin Hosseinzadeh, Saurabh Jha, Anders Jerkstrand, Lindsey Kwok, and Ruediger Pakmor. Many thanks to: Manuel Bautista for providing unpublished collisional strengths for [Ni III] based on Bautista (2001); Lindsey Kwok for sending us the combined spectrum of SN 2021aefx ahead of publication; Elena Sabbi and Brent Tully for confirming the error on the estimated distance to NGC 1566; Ivo Seitenzahl for extracting radial profiles from the 3D delayed-detonation model N100 of Seitenzahl et al. (2013); Steve Shore, the referee of this paper, for his detailed review and valuable comments that led to the detailed comparison of the 3D vs. spherically-averaged yields in Appendix A. This work was supported by the ‘Programme National de Physique Stellaire’ (PNPS) of CNRS/INSU co-funded by CEA and CNES. This research has made use of computing facilities operated by CeSAM data centre at LAM, Marseille, France. D.J.H. acknowledges partial support for this work through the NASA astrophysical theory grant 80NSSC20K0524 and the STScI theory grant HST-AR-16131.001-A. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. This research was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP) which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC-2094 - 390783311. This work made use of the Heidelberg Supernova Model Archive (HESMA; https://hesma.h-its.org) and of the NIST Atomic Spectra Database (https://www.nist.gov/pml/atomic-spectra-database). Note added in proof. Only two days prior to the acceptance of this paper, Kwok et al. (2023b) confirmed the presence of a strong centrally-peaked line due to [Ne II] 12.81 μm in their nebular spectrum of the Type Ia SN 2022pul, which they associate with a violent merger event, in line with the predictions of this paper.

References

Appendix A: 3D vs. spherically-averaged yields

In Table A.1 we compare the total mass and yields of selected isotopes and elements of the original 3D models (M_3D) and their spherically-averaged versions (M_1D) available on HESMA. The spherical averaging causes a systematic overestimate of the total mass by 3–8% and a difference of up to ±10% for several isotopic and elemental yields compared to the original 3D models.

In Table A.2 we repeat this comparison for the original 3D model N100 of Seitenzahl et al. (2013) and 1D models reconstructed from the three orthogonal axes of the 3D Cartesian grid in both positive ({x,y,z}pos) and negative ({x,y,z}neg) directions. We report two sets of quantities, one based on the original (unscaled) density profile, the other based on the density profile rescaled to match the total mass of the original 3D model (1.4 M_⊙). Here the impact on the yields of specific isotopes and elements is in some cases larger than 50%, as is the case for ⁵⁸Ni in the model corresponding to the positive $\widehat{\mathbf{z}}$ direction (+81%), or Ca in the model corresponding to the positive $\widehat{\mathbf{y}}$ direction (+87%).

The yields other than ⁵⁶Ni_{t = 0} are given at 270 d post explosion. In the spherically-averaged models these yields are simply integrated over the ejecta that are used as initial conditions for our 1D radiative-transfer calculations. The 3D yields at 270 d were derived from the final decayed abundances at ∼ 2 Gyr reported in the HESMA abundances.dat files, corrected for radioactive decays based on the abundances of radioactive isotopes at the end of the 3D simulation (∼ 100 s post explosion).

As an example, we consider the iron abundance in the ddt_2013_N100 model, whose final decayed abundance is 7.40 × 10⁻¹ M_⊙. This includes ⁵⁶Fe from the ⁵⁶Ni→⁵⁶Co→⁵⁶Fe decay chain, ⁵⁷Fe from the ⁵⁷Ni→⁵⁷Co→⁵⁷Fe decay chain, and excludes ⁵⁵Fe from the ⁵⁵Co→⁵⁵Fe→⁵⁵Mn. To obtain the abundance at 270 d post explosion, we therefore subtracted the ⁵⁶Co abundance at 270 d (5.81 × 10⁻² M_⊙) as well as the ⁵⁷Co abundance at 270 d (9.48 × 10⁻³ M_⊙), and added the ⁵⁵Fe abundance at 270 d (1.10 × 10⁻² M_⊙), resulting in an iron yield of 6.83 × 10⁻² M_⊙ at 270 d. In principle we would also need to subtract the ⁵⁶Ni and ⁵⁷Ni abundances, but these are essentially zero at 270 d. Likewise, we added the ⁵⁶Co and ⁵⁷Co abundances at 270 d to the final decayed Co abundance to obtain its value at 270 d.

To compute the stable Ni abundance at 270 d, however, we would need to know the initial abundances of several isotopes (such as ⁶⁰Zn and ⁶⁰Cu, which decay to ⁶⁰Ni) that are not included in the HESMA abundances.dat files. Instead, we report the ⁵⁸Ni abundance, which constitutes the dominant stable isotope of Ni (> 60% in the double-detonation models, > 80% in the GCD models, and > 90% in the delayed-detonation and violent merger models; see also Blondin et al. 2022b).

Appendix B: Atomic data

Table B.1 give the number of levels (both super-levels and full levels; see Hillier & Miller 1998 for details) for the model atoms used in the radiative-transfer calculations presented in this paper. N_SL refers to the number of super levels used for the solution of the rate equations, and N_full refers the number of full levels used to solve the transfer equation and compute the observed spectrum. We report the uppermost level for each ion treated in the fourth column. ‘W’ refers to states in which higher ℓ states (usually f or higher) have been combined into a single level. In the last column we give the number of bound-bound transitions in the model ion taking into account all N_full levels. We considered ionisations to and recombinations from the ground state of the next ionisation stage for all elements (i.e. He III, C IV, N IV, O IV, Ne IV, Na II, Mg IV, Al IV, Si V, S V, Ar IV, Ca V, Sc IV, Ti IV, Cr V, Mn IV, Fe VI, Co V, and Ni VI).

Below we give references to the atomic data used (oscillator and collision strengths, with a special attention to forbidden transitions) for all the ions for which we report line identifications in Table 2: Ne II–III, Na I, S II–IV, Ar II–III, Ca II, Ca IV, Fe I–III, Fe V, Co II–IV, and Ni II–IV. We refer the reader to earlier papers for additional details on the atomic data used by CMFGEN (e.g. Hillier 2011, Blondin et al. 2022a).

• Ne II: data for forbidden transitions are from Mendoza (1983);

• Ne III: data for forbidden transitions are very similar (mostly within a few per cent or identical) to those reported in Mendoza (1983);

• Na I: oscillator strengths are from the Opacity Project (Seaton 1987). Collision rates were computed using van Regemorter’s approximation (van Regemorter 1962);

• S II: oscillator and collision strengths for forbidden transitions are from Tayal & Zatsarinny (2010);

• S III: data for forbidden and intercombination (semi-forbidden) transitions are from Tayal (1997), Mendoza (1983), and Huang (1985). Collision strengths are from Tayal & Gupta (1999);

• S IV: data for intercombination (semi-forbidden) transitions are from the compilation of Mendoza (1983). Collision strengths are from Tayal (2000);

• Ar II: ground-state fine-structure cross sections are from Pelan & Berrington (1995). Collision strengths are from Tayal & Henry (1996);

• Ar III: forbidden transition probabilities are from Mendoza & Zeippen (1983). Collision strengths are from Galavis et al. (1995);

• Ca II: unfortunately we were not able to trace the origin of the forbidden transition probabilities for Ca II. The transition probabilities for the ²S_1/2-²D_3/2 (1-2) and ²S_1/2-²D_5/2 (1-3) transitions making up the 7300 Å doublet are both set to A_ul = 1.15 s⁻¹ in our atomic data set, slightly lower than the value reported in the National Institute of Standards and Technology (NIST) Atomic Spectra Database¹¹ (ASD; Kramida et al. 2022), namely A_ul = 1.3 s⁻¹ (Osterbrock 1951). More recent calculations result in significantly lower values (e.g. Kaur et al. 2021 report 0.805 ± 0.007 s⁻¹ and 0.827 ± 0.008 s⁻¹ for the 1-2 and 1-3 transitions, respectively), in line with experimental measurements (0.833 ± 0.007 s⁻¹ and 0.856 ± 0.005 s⁻¹; Barton et al. 2000); we will use these latter experimental values in future work. We ran test calculations using the DDT and MERGER models, which showed a very modest impact. Collision strengths are from Meléndez et al. (2007);

• Ca IV: ground term transition rates are from Nahar & Shafique (2023), as are the collision strengths;

• Fe I: data for forbidden transitions were downloaded from NIST (Kramida et al. 2015) on 7^th September 2016 (the original data are from Brown & Evenson 1995). Collision strengths between the two lowest terms (10 levels) are from Pelan & Berrington (1997);

• Fe II: forbidden transition probabilities among the first 63 levels are from Deb & Hibbert (2011). Collision strengths are from Zhang & Pradhan (1994) and Zhang & Pradhan (1995);

• Fe III: data for forbidden transitions are from Nahar & Pradhan (1996). Collision strengths are from Zhang (1996);

• Fe V: oscillator strengths were computed by R. Kurucz (Kurucz 2009) and obtained through his website¹², complemented with data from NIST¹¹. Transition probabilities for M1 and E2 transitions are from Aggarwal et al. (2017);

• Co II: oscillator strengths were computed by R. Kurucz (Kurucz 2009) and obtained through his website¹² (calculation date: 24^th November 2006). Oscillator strengths for forbidden transitions are from Quinet (1998). Collision strengths are from Storey et al. (2016);

• Co III: oscillator strengths were computed by R. Kurucz (Kurucz 2009) and obtained through his website¹² in 2009. Transition probabilities for forbidden lines in the 3d7 ground configuration are from Hansen et al. (1984). Collision strengths are from Storey & Sochi (2016);

• Co IV: oscillator strengths were computed by R. Kurucz (Kurucz 2009) and obtained through his website¹² in 2009. Transition probabilities for forbidden lines are scaled values for the isoelectronic ion Fe III (Quinet 1996). Collision strengths are also based on Fe III values from Zhang (1996);

• Ni II: oscillator strengths were computed by R. Kurucz (Kurucz 2009) and obtained through his website¹² in 2009. Data for forbidden transitions are from Quinet & Le Dourneuf (1996). Collision strengths are from Bautista (2004);

• Ni III: oscillator strengths were computed by R. Kurucz (Kurucz 2009) and obtained through his website¹² in 2000 (newer calculations are available). Forbidden line transition probabilities are from Garstang (1958). Collision strengths are from this paper (see Sect. 5.1 and Appendix F);

• Ni IV: data for forbidden transitions among low-lying levels are from NIST¹¹ (Kramida et al. 2022; original data from Hansen et al. 1984). Collision strengths are from Fernández-Menchero et al. (2019).

Appendix C: Ionisation fractions

In Figs. C.1–C.7 we show the ratio of the number density of different ionisation states to the total element number density in our reference model set for Ne, S, Ar, Ca, Fe, Co, and Ni.

Fig. C.1. Ratio of the number density of different ionisation states of neon to the total Ne number density in our reference model set.

Fig. C.2. Same as Fig. C.1 but for sulfur.

Fig. C.3. Same as Fig. C.1 but for argon.

Fig. C.4. Same as Fig. C.1 but for calcium.

Fig. C.5. Same as Fig. C.1 but for iron.

Fig. C.6. Same as Fig. C.1 but for cobalt.

Fig. C.7. Same as Fig. C.1 but for nickel.

Appendix D: Statistical estimators

We used several statistical estimators to gauge the fit quality of a given model with SN 2021aefx (see Table D.1). In what follows N_λ is the number of wavelength bins in the observed spectrum, F_obs, i is the observed flux in the i’th wavelength bin, and F_mod, i is the model flux in the same bin (the model spectra were resampled to match the observed wavelength scale).

The median absolute deviation (MAD) is calculated as:

$$\begin{array}{c}\hfill \mathrm{MAD}=\frac{1}{N_{\lambda }}{\displaystyle \sum _{i=1}^{N_{\lambda }}}|F_{\mathrm{mod},i}-F_{\mathrm{obs},i}|.\end{array}$$(D.1)

Since this estimator strongly penalises models that deviate from the observations in a few lines, we can use the logarithm of the flux instead, resulting in:

$$\begin{array}{c}\hfill \mathrm{MADlog}=\frac{1}{N_{\lambda }}{\displaystyle \sum _{i=1}^{N_{\lambda }}}\left|{log}_{10}\left(\frac{F_{\mathrm{mod},i}}{F_{\mathrm{obs},i}}\right)\right|.\end{array}$$(D.2)

The mean fractional error (MFE; see e.g. Vogl et al. 2020) is calculated as:

$$\begin{array}{c}\hfill \mathrm{MFE}=\frac{1}{N_{\lambda }}{\displaystyle \sum _{i=1}^{N_{\lambda }}}\frac{|F_{\mathrm{mod},i}-F_{\mathrm{obs},i}|}{F_{\mathrm{obs},i}}.\end{array}$$(D.3)

Because of the F_obs, i term in the denominator, this estimator gives almost equal weights to regions of low and high flux values. As for the MAD estimator above, we can use the logarithm of the flux to avoid penalising models that deviate from the observations in a few lines:

$$\begin{array}{c}\hfill \mathrm{MFElog}=\frac{1}{N_{\lambda }}{\displaystyle \sum _{i=1}^{N_{\lambda }}}|{log}_{10}\left(\frac{F_{\mathrm{mod},i}}{F_{\mathrm{obs},i}}\right)/{log}_{10}(F_{\mathrm{obs},i})|.\end{array}$$(D.4)

The ‘Score’ estimator is defined in Omand & Jerkstrand (2023). Here we used its mean value:

$$\begin{array}{c}\hfill \mathrm{Score}=\frac{1}{N_{\lambda }}{\displaystyle \sum _{i=1}^{N_{\lambda }}}{\left[{log}_{10}(\frac{F_{\mathrm{mod},i}}{F_{\mathrm{obs},i}}.)\right]}^2\end{array}$$(D.5)

The logarithm ensures that models that deviate from the observations in a few lines are not overly penalised, but more so than in the previous MADlog and MFElog estimators due to the squared exponent.

Appendix E: Relativistic effects

The prominent [Ar III] 8.99 μm line displays a blue-to-red tilt in its flat-top profile. As noted previously, our ‘flat’-top line profiles are tilted the opposite way due to relativistic beaming (Fig. E.1), which results in a 2D ≈ 3% contrast between the blue and red edges at ±4500 km s⁻¹ from the rest wavelength, where D is the so-called Doppler factor:

$$\begin{array}{c}\hfill D=\frac{1}{\gamma (1-\beta cos\theta )},\end{array}$$(E.1)

Fig. E.1. [Ar III] 8.99 μm line profile in our DDT model with (blue) and without (red) relativistic corrections in the observer-frame calculation of the spectrum. The inset shows a close-up view of the ‘flat’-top portion of the profile, highlighting the reversal of the tilt direction.

where β = v/c and $\gamma =1/\sqrt{(1-{\beta }^2)}$ is the Lorentz factor. While line overlap can conspire to reverse the direction of this tilt in some cases, this does not seem to apply to the [Ar III] line. The observed ∼ 10% increase in flux from blue to red then suggests a pronounced asymmetry in the argon distribution (see discussion in DerKacy et al. 2023).

Appendix F: Collision strengths for Ni III

We report transition probabilities (A_ul) and effective collision strengths (Υ_lu) in the temperature range 1000–40 000 K among the lowest ten levels of Ni III (see Table F.1) in Table F.2.

Axelrod (1980) proposed general approximations for collision strengths for forbidden-line transitions based on [Fe III], which depend on the product of the statistical weights of the lower and upper levels (g_lg_u):

$$\begin{array}{c}\hfill {\mathrm{{\rm Y}}}_{\mathit{lu}}=a_{\mathrm{IR}}w\frac{g_l{g}_u}{8},\end{array}$$(F.1)

where w ≈ 0.03 and a_IR = 6 for infrared transitions (≳10 μm), and = 1 otherwise. As shown in Fig. F.1, these relations only provide a very crude estimate for most transitions, including among low-lying levels.

Fig. F.1. Effective collision strength (Υlu) at 5000 K for forbidden transitions of Ni III versus the product of the statistical weights of the lower and upper levels. Different symbols are used for transitions corresponding to wavelengths < 10 μm (filled circles) or ≥10 μm (open diamonds). We also highlight transitions among the first ten levels (see Table F.1). Overplotted are the approximations of Axelrod (1980) for the two wavelength regimes. The Υlu values deviate by more than one order of magnitude for many transitions (including for some among low-lying levels).

Appendix G: Spectra for the complete gcd and doubledet model sets

We show the synthetic spectra for the complete set of pulsationally assisted gravitationally-confined detonation models and double-detonation models compared to SN 2021aefx in the range 0.35–14 μm in Figs. G.1 and G.2, respectively.

Fig. G.1. Similar to Fig. 4 but for the full set of pulsationally assisted gravitationally-confined detonation models. The gcd_2021_r10_d1.0 (GCD) model is used as the reference for this class in this paper.

Fig. G.2. Similar to Fig. 4 but for the full set of double-detonation models. The doubledet_2021_M1002_1 (DBLEDET) model is used as the reference for this class in this paper.

All Tables

All Figures

	Fig. 1. Density profiles at 270 d post explosion for our reference model set: the delayed-detonation model ddt_2013_N100 (DDT), the pulsationally assisted gravitationally-confined detonation model gcd_2021_r10_d1.0 (GCD), the double detonation model doubledet_2021_M1002_1 (DBLEDET), and the violent merger model merger_2012_11+09 (MERGER).
In the text

	Fig. 2. Abundance profiles at 270 d post explosion for our reference model set.
In the text

	Fig. 3. Characteristic ejecta properties at 270 d for our reference model set. Top: temperature profiles. Middle: electron density profiles. Bottom: mean Fe ionisation profiles. The mean ionisation is defined as ∑iini+/∑ini+, where ni+ is the number density of ionisation stage i for Fe, such that a mean ionisation of ∼2 indicates that Fe2+ is the dominant stage.
In the text

Fig. 4. Spectra of our reference model set at 270 days post explosion compared to SN 2021aefx over the wavelength range 0.35–14 μm (logarithmic scale), with flux in Fν on a non-linear (arcsinh) scale, as in Fig. 1 of Kwok et al. (2023a). The SN 2021aefx spectrum has been corrected for redshift and extinction (assuming a host-galaxy reddening of 0.097 mag and a MW reddening of 0.008 mag). The synthetic fluxes correspond to the same assumed distance to SN 2021aefx of 18 Mpc; they have not been rescaled or normalised in any way. We include selected line identifications based on their maximum Sobolev EW (see Table 2). The insets show normalised line profiles in velocity space for selected transitions, illustrating variations in line widths and morphology.

In the text

	Fig. 5. Normalised cumulative integral of the flux per unit wavelength over the range 0.35–14 μm in SN 2021aefx and in our reference model set. The vertical lines mark the boundaries between the optical-NIR and NIR-MIR ranges as defined in Sect. 3.3.
In the text

Fig. 6. Same as Fig. 4 but for restricted wavelength ranges (from top to bottom): optical (0.35–1 μm), NIR (1–5 μm), lower MIR (5–14 μm) and upper MIR (14–28 μm). The ℱX−Y μm label gives the fraction of the total optical to MIR flux (0.35–14 μm) for SN 2021aefx emitted within the wavelength range of each plot. We include line identifications based on their Sobolev EW (see Table 2). Transitions connected to the ground state are marked with a ‘†’ symbol. For consecutive lines within a multiplet (connected by a horizontal line), we give the gf-weighted mean wavelength of the transitions. Lines that only appear in one model class are labelled with the corresponding colour.

In the text

	Fig. 7. Density profiles for the spherically-averaged ddt_2013_N100 (DDT) model and along the three orthogonal axes of the original 3D Cartesian grid, in both positive ({x,y,z}pos) and negative ({x,y,z}neg) directions. The bottom panel shows the density ratio with respect to the DDT model.
In the text

	Fig. 8. Abundance profiles for Ni, Co, Fe, Ca, Ar, and S in the spherically-averaged ddt_2013_N100 (DDT) model and along the three orthogonal axes of the original 3D Cartesian grid, in both positive ({x,y,z}pos) and negative ({x,y,z}neg) directions. The lower panels show the logarithm of the abundance ratio with respect to the DDT model.
In the text

	Fig. 9. Similar to Fig. 4 but for the ddt_2013_N100 (DDT) model along the three orthogonal axes of the original 3D Cartesian grid, in both positive ({x,y,z}pos) and negative ({x,y,z}neg) directions. We also show the original (spherically-averaged) DDT model.
In the text

Fig. 10. Impact of different values of the effective collision strengths amongst the lowest 3F term in Ni III (see Table 3) on the lines at 7.35 μm and 11.00 μm in our reference model set. In the top row we show the DDT and MERGER models where these lines are largely overestimated. In the bottom row we show the GCD and DBLEDET models where the predicted strengths for both lines are comparable to those seen in SN 2021aefx (note the smaller ordinate range compared to the top row).

In the text

	Fig. 11. Impact of the minimum neon mass fraction on the [Ne II] 12.81 μm line in the GCD model. We can reproduce the observed 12.8 μm feature in SN 2021aefx by artificially setting the minimum Ne mass fraction to 10−2 throughout the ejecta.
In the text

Fig. 12. Similar to Fig. 4 but for the ddt_2013_N100 (DDT) model with different values for the volume-filling factor f used to approximate a uniformly-clumped ejecta (see Wilk et al. 2020). We show the original (unclumped) DDT model and clumped models with f = 0.75, 0.5, 0.25, 0.1, and 0.05. We highlight lines of neutral ions that emerge in the most clumped models, in particular a strong [Ni I] 3.12 μm line.

In the text

	Fig. 13. Total integrated luminosity in all lines of Ni I–IV in the wavelength range 0.35–14 μm versus stable Ni mass in our complete model set. There is a strong correlation between both quantities, with one notable outlier (the MERGER model). The reference delayed-detonation model (DDT) is highlighted with a grey circle.
In the text

	Fig. 14. Left: [Ar III] 8.99 μm line profiles in the MCh DDT model compared to sub-MCh double-detonation models with varying He-shell masses (the reference DBLEDET model has the lowest He-shell mass of 0.02 M⊙). Right: Ar distribution, highlighting the larger extent of the Ar hole in these sub-MCh models compared to the MCh DDT model.
In the text

	Fig. C.1. Ratio of the number density of different ionisation states of neon to the total Ne number density in our reference model set.
In the text

	Fig. C.2. Same as Fig. C.1 but for sulfur.
In the text

	Fig. C.3. Same as Fig. C.1 but for argon.
In the text

	Fig. C.4. Same as Fig. C.1 but for calcium.
In the text

	Fig. C.5. Same as Fig. C.1 but for iron.
In the text

	Fig. C.6. Same as Fig. C.1 but for cobalt.
In the text

	Fig. C.7. Same as Fig. C.1 but for nickel.
In the text

	Fig. E.1. [Ar III] 8.99 μm line profile in our DDT model with (blue) and without (red) relativistic corrections in the observer-frame calculation of the spectrum. The inset shows a close-up view of the ‘flat’-top portion of the profile, highlighting the reversal of the tilt direction.
In the text

Fig. F.1. Effective collision strength (Υlu) at 5000 K for forbidden transitions of Ni III versus the product of the statistical weights of the lower and upper levels. Different symbols are used for transitions corresponding to wavelengths < 10 μm (filled circles) or ≥10 μm (open diamonds). We also highlight transitions among the first ten levels (see Table F.1). Overplotted are the approximations of Axelrod (1980) for the two wavelength regimes. The Υlu values deviate by more than one order of magnitude for many transitions (including for some among low-lying levels).

In the text

	Fig. G.1. Similar to Fig. 4 but for the full set of pulsationally assisted gravitationally-confined detonation models. The gcd_2021_r10_d1.0 (GCD) model is used as the reference for this class in this paper.
In the text

	Fig. G.2. Similar to Fig. 4 but for the full set of double-detonation models. The doubledet_2021_M1002_1 (DBLEDET) model is used as the reference for this class in this paper.
In the text