Type Ia supernovae in the age of JWST: Finding the ‘right’ questions and the path to answers

Understanding the Physics of thermonuclear explosions of a White Dwarf star (WD), so called Type Ia Supernovae (SNe Ia), provide a playground for modern physics, computational methods, and are a key to modern cosmology. We identify new and investigate a variety of observational signatures of underlying physical processes related to the thermonuclear runaway, the flame propagation and the environment. Being intrinsically multi-dimensional phenomena, probing the physics requires multi-dimensional radiation-hydrodynamics and MHD simulations. For this task, we developed and employed methods for photon transport for the X-, gamma- and of low energy and of positrons under non-LTE conditions. We identify signatures in the light curves and spectra, in particular, line profiles and polarization spectra. Consistent treatment of high energy processes is critical. Therefore, our framework and results can be used directly a variety of scenarios for SNe Ia including merging WDs and explosions of sub-Chandrasekhar mass WDs. Current simulations have limitations but, nevertheless, when combined with recent JWST and VLT observations solutions emerge to many of decade old problems on the ignition process, flame physics and thermonuclear explosion.


Introduction
Type Ia supernovae (SNe Ia) explosively destroy C/O white dwarfs (WDs) in binary or triple systems and produce 50 % of the iron-group elements in the Universe.They have gained wide attention because the luminosity decline rate relation led to the discovery of Dark Energy.But understanding the directional dependence of the luminosity is needed for precision cosmology and to test the underlying fundamental physics.The brilliance of the top-level discoveries outshines the enormous diversity of SNe Ia that has emerged in recent years and led to an even larger diversity of three main explosion scenarios [1]: a) One of the leading explosion scenario considers a surface He detonation (HeD) that triggers a detonation of a sub-M Ch with a low-density core C/O core [2,3,1].b) A second scenario considers delayed detonation (DD) [4], in which the explosion of a near-M Ch -mass WD is triggered by compressional heating in high densities resulting in the production of electron-capture elements (EC).Here, the flame propagation starts as a deflagration and transitions to a detonation.c) In the third leading scenario, two WDs merge or collide, possibly head-on, in a triple system [5], causing large-scale asymmetries and low-density burning.Despite all efforts, none of the current modeling efforts is 'first-principle', and most simulations of light curves and spectra assume spherical geometry, often on angle averaged quantities.For a collection of reviews by all groups, see [6].With tuning, all models can reproduce overall light curves and total-flux spectra because they depend on the same network of nuclear reactions and similar WD structures, abundances, and explosion energies, a fact described as "stellar amnesia".By contrast, geometry encapsulates explosionand progenitor-specific information (Fig. 1).Geometry is sensitive to the companion star and the  [7].Middle: Temperature and turbulent velocity structure at the onset of thermonuclear runaway of the inner 100 km of a WD [8].Right: Helium-triggert detonation in a sub-MCh WD [9].Lower left: Developing RT-instabilities in deflagration fronts result in typical scales of ∼ 1000 km/sec during the homologous expansion phase [10].Middle: Influence of the magnetic field on the nuclear burning front in a flux tube of 240 km at 0.6 sec for WD-typical conditions.We show the burning rate in arbitrary units in 6 panels for B between 0 and 10 12 G (bottom to top); all other parameters are identical.B-fields ≤ 10 6 G are inferred from late time light curves [11].Right: 56 N i abundance as a result of an off-center deflagration-detonation transition (black dot) [12].circumstellar environment(CSM), the physics of the ignition and flame propagation, and a new generation of models make clear and new qualitative and quantitative predictions.Moreover, features of EC-elements are highly blended in the optical and near IR, making JWST with MIR a game-changer.For the first time, it is possible to obtain the full 3-D image of SNe Ia by combining multi-messenger and time-domain observations from days to years after the explosion: optical to NIR total-flux and polarization spectra probe C/O, O/Ne/Mg, and Si/S, nebular optical/NIR spectra cover Si/S/Ca/Co/Fe, and mid-IR spectra Mn and Cr/Ar/Fe/Co/Ni.However, the interpretation requires a new generation of multi-dimensional non-LTE radiationhydro simulations which, still are challenging and limited by resolution as described below.

Numerical Methods
Simulations for transients requires finding the simultaneous solutions of equations for nuclear networks, the statistical equations needed to determine the atomic level population, treat the flame physics, the equations of state, the opacities, and the hydro and radiation problems.Our HYDrodynamical RAdiation HYDRA code consists of physics-based modules to provide a solution (Fig. 2).The basic set of equations for solving the individual physical problems can be found in [13,14], and references therein.Consistency between the solutions is achieved iteratively [15,16].The second approach employed is the extended use of hybrid-schemes of grid-based and Monte-Carlo methods with detailed physics problems solved by MC.Hydrodynamics: The hydro modules use an explicit Piecewise Parabolic Method (PPM) by [21] to solve the compressible reactive flow equations with variable adiabatic gradients.PPM is implemented as a step followed by separate remaps of the thermal and kinetic energy to avoid numerical generation of spurious pressure disturbances.Subsonic reactive fronts are treated [22], based on 3D hydrodynamical simulation.For 3D geometries, [10] in case of spherical symmetry or, in 3D,the burning operator is used as given in equations 8 & 9 of [23].High-density equation of state and nuclear reactions: The EOS is for a partially degenerate and partially relativistic Bose-and Fermi-gas plus interactions due to effects Coulomb corrections, quantum-relativistic effects on the electron component, and electron-positron pair production and electron screening.For nuclear reactions, we use an extended network of 218 isotopes based on the implementation by [24] but with modified matrix solvers and updated, parameterized reaction rates , and employ nuclear statistical equilibrium for temperatures in excess of 6.5×10 9 g/cm 2 .All relevant strong, weak and photon interaction have been taken .The   [15].The spectra are based on our spherical and full 3-D radiation transport modules which 910 (67/67/67) depth points, 20,000 (2,000) frequency points, and 520 (50) non-LTE-levels, respectively.Even today, many simulations are spherical models with low R (∼50 to 100) [17].The differences were ≈ 20% in flux at 4000 Å and B − V color, which is comparable to the directional dependence of SN2004dt derived with a spatial grid of 300 2 points ( [18], Fig. 4).Multi-D needs high resolution explaining why multi-D non-LTE comes to age only now.Middle: The effect of the spatial R on the mean intensity relative to a black body is shown using n(r) depth points on an equidistant mesh for conditions resembling an SN Ia density structure at maximum light, namely a scattering optical depth of τ = 30 with a thermalization fraction of 10 −3 [19,16].J/B = 1 is required to force LTE-population numbers, and errors of a factor of 3 are found for for n(r)=100 while n(r) ≈ 300 is actually required [12].Right: As part of an ongoing PhD project for the VET solver, the domain deposition is currently being implemented by Mera Evans using 'pancake slicing' in the spatial domain and cycling between the three Cartesian coordinates.Domain splitting and initialization are shown before bottom-to-top directional sweep [20].
individual rates between isotopes are fitted by expansions in ρ and T with factors tabulated in the reaction-network library REACLIB by Thielemann, and modified weak reactions published by [25].The expansion of the rates assume a Boltzman distribution for the energies of the particles.Alternatively, the rates are calculated by integration of the cross-sections using a partial degenerate Fermi-gas.Low density equation of state, opacities and source functions: For the atomic models, we reduce the number of energy levels by use of level-merging/super-levels, with the assumption that coupling between the merged levels being in thermodynamical equilibrium.This implies full frequency redistribution of individual transitions (see below).For the radiation transport.For Monte Carlo transport, in the simulations shown here, we use narrow line limit [26] and take into account the probability for absorption along the way by lines of higher frequency.The module for molecular kinetics use rates [27].For the atomic level population, we solve the full set of statistical equations, namely for the rate, charge and particle conservation.The time-dependence in the rate equations We solve the full set of statistical equations including particle and charge conservation to determine the population.The non-thermal excitation by hard γ-rays and electrons are added into the radiative rates R and collisions with fractions according to the individual level densities following the cascading in energy by MC [28,12].Time scales in the rate equations are dominated by the slow bound-free transitions.Therefore, time-dependent solution for the ionization balance can be obtained by solving an inhomogeneous ODE analytically or as a simple system of linear equations [29].Coupling of radiation transport, statistical and hydro equations: We employ to acceleration schemes: a) Accelerated lambda iteration (ALI, e.g.[30]); b) equivalent-2-level approach [31,32] which drive non-locality in the spatial and in the frequency space via levels to capture both the high-optical depth individual lines and global frequency space needed, e.g. in nebular spectra [23].Radiation transport for low energy photons: For the time-dependent transport, we use variable Eddington Tensor solvers with full Monte-Carlo simulation in space and frequency [33,15,34,16].For the 3D case, the variable Eddington-Tensor method is used which is accurate to the order O(u/v) and we neglect acceleration terms [35,36,33,16,37].The system of equations is given by with the tensor T and, following Auer (1971), the variable q ν as defined by: and ) Ω dΩ, and P ν ( x) = 1/c I ν ( x, Ω) Ω ⊗ Ωd Ω.
(5) with the integrals over ν denote the corresponding radiative components in the hydro-equations.For 3D, the advection frequency redistribution terms enter via the MC solution.For some tests of 3D vs. spherical solutions (Fig. 2).Again, we use stationary transport for the closure relation and tested the result against the spherical transport module.Monte Carlo (MC) methods are used including advection and aberration for the closure relation of low energy photons for flux and polarization spectra, and similar to hard γ and positron transport [38,39,40,41], and the cascading of high-energy thermal leptons [12].Beyond velocity gradients, the redistribution of energy in frequency space is caused by a large number of overlapping lines are important.E.g., in the optical, a typical of 10 6 transitions are taking into account over ≈ 10 4 Å.This leads to a coupling of photons in the frequency space via the rate-equations.The same MC transport used for the closure allows to calculate the redistribution functions ψ ν in the source terms numerically by weighting the neighboring frequencies with the overlap of a specific transition so that the frequency integral over S ν is conserved [29].

Some applications and first results based on JWST and SPECPOL
Up to now, all 5 'classical' normal thermonuclear SNe with MIR spectra show nebular spectra dominated by forbidden iron-group elements and multiple stable 58 N i, indicative for high-density burning (Fig. 1) found in M Ch [50,51,52,46] including one underluminous SN2022xkq (Fig. 5).The small Doppler shift of N i suggests close to center, single spot ignition.Pre-existing turbulent velocities are consistent with the line-width [53,48] whereas off-center, multi-spot off-center ignitions would produce extended mixing by RT-instabilities and, thus, broad lines [49].The lack of microscopic mixing of Ni possibly due to high B-fields [11] is evident in SN2022aefx and SN2022xkq based on our simulations.Line-asymmetries in Ar and Co are consistent with off-center delayed-detonation transitions in the regime of distributed burning from both spectropolarimetry (Fig. 4) and line profiles [23,54].Our current simulations have resolutions of R ≈ 300 posing a severe restriction for probing small scales.Detailed results of the JWST MIR-spectra and spectropolrimetry programs can be found in [12,52,54,46].[45,46] obtained with MIRI/MRS at JWST [47] and models with ρc = 2 × 10 9 g cm −3 (solid red) and 5 × 10 8 g cm −3 (dotted blue) including line identification based on the simulations.All spectra have been normalized using a distance module of 32.52 mag.The JWST pipeline provides excellent calibration and background reduction, a combination of observations and synthetic spectra is needed to properly reduce the background if it dominates the flux as in Channel 4 of JWST.The correction log of the calibration factor is given in cyan.The simulations allow to analyze the data on a 10% level.The low central density model strongly over-predicts the flux in the 7 and 9 µm [Ar] providing tight limit on the central density of the exploding WD and the presence of narrow N i lines, somewhat surprisingly, identify the as a MCh type rather than a sub-MCh, helium triggered explosion.The narrow line strongly suggest single-point and mostly central ignition [8,48] rather than a strong off-center multi-spot ignition [49].

Figure 1 :
Figure1: Progenitor structure and multi-dimensional imprints.Upper left: Central WD density as a function of mass and the influence of element production[7].Middle: Temperature and turbulent velocity structure at the onset of thermonuclear runaway of the inner 100 km of a WD[8].Right: Helium-triggert detonation in a sub-MCh WD[9].Lower left: Developing RT-instabilities in deflagration fronts result in typical scales of ∼ 1000 km/sec during the homologous expansion phase[10].Middle: Influence of the magnetic field on the nuclear burning front in a flux tube of 240 km at 0.6 sec for WD-typical conditions.We show the burning rate in arbitrary units in 6 panels for B between 0 and 10 12 G (bottom to top); all other parameters are identical.B-fields ≤ 10 6 G are inferred from late time light curves[11].Right: 56 N i abundance as a result of an off-center deflagration-detonation transition (black dot)[12].

Figure 2 :
Figure 2: Block diagram of our numerical scheme to solve radiation hydrodynamical problem.

Figure 3 :
Figure3: The need for high-resolution multi-D simulations.Left: Influence of the discretization on frequency, dimensions, and the complexity of the atomic models based on a spherical sub-luminous DDT model as an example to demonstrate the importance of resolution (R) and complex non-LTE models, reflecting our status of 2002[15].The spectra are based on our spherical and full 3-D radiation transport modules which 910 (67/67/67) depth points, 20,000 (2,000) frequency points, and 520 (50) non-LTE-levels, respectively.Even today, many simulations are spherical models with low R (∼50 to 100)[17].The differences were ≈ 20% in flux at 4000 Å and B − V color, which is comparable to the directional dependence of SN2004dt derived with a spatial grid of 300 2 points ([18], Fig.4).Multi-D needs high resolution explaining why multi-D non-LTE comes to age only now.Middle: The effect of the spatial R on the mean intensity relative to a black body is shown using n(r) depth points on an equidistant mesh for conditions resembling an SN Ia density structure at maximum light, namely a scattering optical depth of τ = 30 with a thermalization fraction of 10 −3[19,16].J/B = 1 is required to force LTE-population numbers, and errors of a factor of 3 are found for for n(r)=100 while n(r) ≈ 300 is actually required[12].Right: As part of an ongoing PhD project for the VET solver, the domain deposition is currently being implemented by Mera Evans using 'pancake slicing' in the spatial domain and cycling between the three Cartesian coordinates.Domain splitting and initialization are shown before bottom-to-top directional sweep[20].

Figure 4 :
Figure 4: An off-center DDT model for a normal-bright SNe Ia viewed at an inclination of +45 • compared to SN2019np observed by our SPECPOL VLT/FORS program on days −14.5, −11.4,−6.4,+5 and +14.5 w.r.t.maximum light [42].Left: the evolution of the polarization spectra.Right, upper: Continuum polarization at days −14.5 (middle) and +14.5 (right) with the colored areas indicating the value consistent with the observations.Right, middle: Total-flux spectrum and photospheric radius at maximum.The observed p is reproduced both across SiII λ6360 and in the continuum.Many of the weak features in the polarization spectra (left) can be attributed to line transitions.SN 2019np was an overall spherical explosion with an asymmetric surface layer of ∼ 5 10 −3 M and an aspherical 56 Ni core.Polarization and total-flux spectra form in the same part of thev extended atmosphere (right middle vs. insert in upper middle panel).A comparison between an off-center DDT model for SN2004du (lower right) shows similar features but much higher p [43] in lines though, still, a low continuum p. Within this framework, this can be understood such that in SN 2004du the DDT was triggered at MDDT = 0.25M (vs.MDDT = 0.5M in SN 2019np) and at an optimal viewing angle, namely +30 • .Lower, middle: The unblended nebular [Fe II] line profile at 1.644 µm on day +210 viewed from −90 o with MDDT = 0., 0.3 and 0.9 MWD in comparison to the underluminous SN2020qxp (black) observed at Keck (PI: C. Ashall) shows that nebular spectra offer complementary means to determine MDDT and the viewing angle[44].Discrepancies in the results would hint towards problems with the explosion scenario employed.

Figure 5 :
Figure 5: Comparison between the MIR spectra of the subluminous SN2022xkq[45,46] obtained with MIRI/MRS at JWST[47] and models with ρc = 2 × 10 9 g cm −3 (solid red) and 5 × 10 8 g cm −3 (dotted blue) including line identification based on the simulations.All spectra have been normalized using a distance module of 32.52 mag.The JWST pipeline provides excellent calibration and background reduction, a combination of observations and synthetic spectra is needed to properly reduce the background if it dominates the flux as in Channel 4 of JWST.The correction log of the calibration factor is given in cyan.The simulations allow to analyze the data on a 10% level.The low central density model strongly over-predicts the flux in the 7 and 9 µm [Ar] providing tight limit on the central density of the exploding WD and the presence of narrow N i lines, somewhat surprisingly, identify the as a MCh type rather than a sub-MCh, helium triggered explosion.The narrow line strongly suggest single-point and mostly central ignition[8,48] rather than a strong off-center multi-spot ignition[49].