Spectral synthesis techniques for supernovae and kilonovae

TL;DR

This paper reviews techniques for modeling the spectra of supernovae and kilonovae, focusing on computational methods and their evolution.

Contribution

The paper provides a comprehensive review of spectral synthesis techniques and their computational challenges for supernovae and kilonovae.

Findings

Spectral synthesis modeling has evolved from stellar winds to supernovae and kilonovae.

Current codes use various approximations for central physical processes.

Similarities and differences in numeric schemes are identified for improved models.

Abstract

Supernovae (SNe) and kilonovae (KNe) are the most violent explosions in cosmos, signalling the destruction of a massive star (core-collapse SN), a white dwarf (thermonuclear SN) and a neutron star (KN), respectively. The ejected debris in these explosions is believed to be the main cosmic source of most elements in the periodic table. However, decoding the spectra of these transients is a challenging task requiring sophisticated spectral synthesis modelling. Here, the techniques for such modelling is reviewed, with particular focus on the computational aspects. We build from a historical review of how methodologies evolved from modelling of stellar winds, to supernovae, to kilonovae, studying various approximations in use for the central physical processes. Similarities and differences in the numeric schemes employed by current codes are discussed, and the path towards improved models is laid out.

Figures

Characteristic evolution of SNe and KNe over the three phases of diffusion phase, early tail phase, and late tail phase

Full text

Introduction

Historical background and overview

Explosive transients have always been one of the cornerstones of astronomy. Galactic supernovae have brought attention to the otherwise quiescent night sky for as long as human civilization has existed, with the oldest surviving records being those of Chinese astronomers observing SN 185 over 1800 years ago.

It was the Swedish astronomer Knut Lundmark, who in his 1925 treatise “The motions and distances of spiral nebulae” was the first to realize that there was a particular subgroup of novae that appeared fundamentally different to the usual ones, being much brighter. Lundmark refered to these as giant novae. The name that would come to stick, however, was supernovae, first used in Baade and Zwicky (1934). Baade and Zwicky speculated that supernovae represented the collapse of stellar cores to neutron stars, a remarkably apt inferrence coming just two years after the discovery of the neutron by James Chadwick in 1932. Ironically, all the supernovae that were known at that time were of the Type Ia variant, representing a completely different phenomenon; the thermonuclear destruction of a white dwarf. Today, we know that the collapse of the cores of massive stars to neutron stars (“core-collapse supernovae”) make up about 80% of all supernovae by unit volume, with the Type Ia class the other 20% (Li et al. 2011). But it was not until the mid-1980s that the various observational classes could be correctly associated with the right type of explosion (Wheeler and Levreault 1985; Filippenko and Sargent 1986; Gaskell et al. 1986).

The key to that development, and much of supernova science since, was spectral observations and modelling of these spectra. One of the first papers (a conference proceeding) which presented supernova synthetic spectra was Branch (1980). He was the first to develop and adapt the Schuster–Schwarzschild modelling approach (a frequency-independent photophere with a scattering atmosphere outside) to supernovae, identifying the use of important approximations such as the Sobolev formalism for line transfer. With this tool in place, quite detailed comparisons between model spectra of white dwarf explosion simulations and observations of Type I SNe (Fig. 1) could be carried out starting with the works of Branch et al. (1983, 1985). The tool built by David Branch and his group would eventually become the famous SYNOW code, which is available at https://c3.lbl.gov/es/. It is still frequently used today for line identifications and rapid model investigations.Fig. 1. Example of an early SN spectral model (bottom) for the photospheric spectrum of a white dwarf explosion model compared to an observed spectrum (top). From this modelling, the first identification of which lines are important for the different observed featured could be established. Image reproduced with permission from Branch et al. (1985), copyright by AAS

As the ejecta expand the photosphere eventually disappears and the inner ejecta become visible, emitting nebular emission lines. The first step for computing spectral models in this phase was taken also in 1980, in the remarkable PhD thesis by Timothy Axelrod at Berkeley, supervised by Tom Weaver and Stan Woosley (Axelrod 1980). This work lays the foundation for the various pieces of physics going into such models, including non-thermal and NLTE (Non-Local Thermodynamic Equilibrium) physics. It still today sets the standard for many aspects of nebular-phase spectral modelling. The theoretical foundation by Axelrod has been used in several later nebular-phase codes, e.g. Mazzali et al. (2001); Maeda et al. (2006).

The explosion of SN 1987A brought more workers into the field. In Stockholm, Claes Fransson developed spectral models for its late emission (Fransson and Chevalier 1987), and then extended this to Type Ib supernovae (Fransson and Chevalier 1989). In London, Leon Lucy developed a Monte Carlo code for Schuster–Schwarschild modelling (Lucy 1987). The Lucy code was further developed and applied in Ruiz-Lapuente et al. (1992), Mazzali and Lucy (1993), Mazzali (2000).

Starting around 2005, significant developments took place for supernova spectral synthesis techniques, much of it inspired by the works of Leon Lucy over the preceding years to develop Monte Carlo techniques beyond the first simplified uses (Lucy 1999, 2003, 2005). The 3D LTE Monte Carlo codes SEDONA (Kasen et al. 2006) and ARTIS (Kromer and Sim 2009) appeared. The 1D NLTE codes (with radiative transfer) SUMO (Jerkstrand et al. 2011, 2012), NERO (Maurer et al. 2011) and CMFGEN (Hillier and Dessart 2012, describe its adaptations to SN applications) were developed. A first public Schuster–Schwarzschild code, TARDIS, was released (Kerzendorf and Sim 2014).

Just as much of supernova modelling became an extension of methods, tools and concepts originally devised for stellar winds and H I regions, mostly during the 1970s (e.g. Lucy and Solomon 1970; Dalgarno and McCray 1972), it in turn became the spring-board for kilonova modelling. This era started in earnest with the paper of Metzger et al. (2010), where the supernova code SEDONA was extended to be used for kilonovae. Other 3D LTE Monte Carlo codes capable of KN modelling were developed by Tanaka and Hotokezaka (2013, the code is not formally named, we will refer to it as TH13here) and Wollaeger et al. (2013, SuperNu). In this last paper, new technical steps were taken by the development and application of implicit Monte Carlo methods. The code POSSIS (Bulla 2019, 2023) operates on similar principles as ARTIS and TH13, and the SN code JEKYLL (Ergon et al. 2018) implements several computational efficiency improvements to the Lucy method. The SUMO code was recently adapted to KN modelling (Pognan et al. 2022b), which has opened up the path towards NLTE modelling of KNe. The first steps towards 3D NLTE modelling, so far for supernovae only, have also been taken (Botyánszki et al. 2018; Shingles et al. 2020; van Baal et al. 2023).

Looking back at the tapestry of developments, one can see how astrophysics grows by small extension steps into adjacent, related areas. That process takes time, as in decades. It is interesting to contextualize that “application diffusion” process against changes in computing power. Sometimes growth in computing power can drive such steps to be taken. But more often, they tend to be application-driven. For the latter case, it may happen that methods used eventually get outdated, with the original derivations occurring under the constraints of computing power limitations orders of magnitudes stricter than today. It is worth contemplating these aspects whenever a code or methodology is studied; sometimes simplifications and approximations done in the original model are no longer necessary, and modern computing power can be taken advantage of to improve upon the approach. This is one of several things we gain when studying the history and gradual development of methodologies.

Classification

Supernovae come in two main types; the explosion of massive stars, called core-collapse supernovae, and the explosion of white dwarfs, called thermonuclear supernovae. The resulting SNe are classified in a system using both spectral and light curve properties (Filippenko 1997). Presence of H lines in spectra leads to a “Type II” classification, and absence to a “Type I” classification. For the Type I class, presence of He lines leads to “Type Ib” subclass, and absence to “Type Ic”. A particular group of Type Ic SNe have unusually broad lines, and are referred to as “Type Ic-BL” SNe. The Type II’s are divided into subclasses of Type IIP if the light curve has a plateau, and Type IIL if it instead has a linear decline. A rise and decline, like SN 1987A, gives a Type II-pec labelling. Finally, if narrow H lines are seen (thought to arise from the circumstellar medium rather than the SN), the classification is Type IIn.

The thermonuclear supernovae are classified as “Type Ia”. The destruction of the white dwarf involves ignition of its large content of carbon. Because conditions are degenerate, such nuclear burning meets no damping by pressure expansion, and a runway occurs. The ignition can occur either when the white dwarf accretes matter from a companion star, or when it merges with another white dwarf. The smaller variation possible both for the exploding star (a white dwarf close to the Chandrasekhar mass), and the resulting nucleosynthesis (everything is burnt to the iron-group), means Type Ia SNe show more uniform properties and correlations than CCSNe. This allowed them to be the standard candle tool (Phillips 1993) that eventually led to the discovery of the accelerating expansion of the Universe (Riess et al. 1998; Perlmutter et al. 1999). However, rare subclasses also exist (Taubenberger 2017).

While SNe have been observed and scientifically studied for over a century, kilonovae were discovered only in 2017. These transients involve the merger of two neutron stars, or a neutron star and a black hole—a much more rare phenomenon than the death of a massive star or a white dwarf. While most of the mass in the merging bodies will form a black hole, about 1% or so ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 0.05\, M_\odot $$\end{document} ) will be ejected from the gravitation potential. This material has gone through the r-process—rapid capture of neutrons forming heavier and heavier elements. The composition is therefore trans-iron elements. Many of the isotopes formed are radioactive, which provides an important power source to generate a bright transient from the ejected material. Because the mass is low, but the velocity high (about 10% the speed of light), KNe rise and decline over just a few days, compared to weeks or months for SNe. With only one solid detection, and a handful of candidates, no classification system yet exists for them.

Motivation and goals of supernova and kilonova spectral synthesis

Characteristic evolution of SNe and KNe over the three phases of diffusion phase, early tail phase, and late tail phase

Figure 2 shows a schematic illustration of the evolution of explosive transients, and what layers spectra probe at different phases. As optical depths decline with time, it becomes possible to see deeper and deeper into the nebula, and therefore is time-sequencing important to obtain full-ejecta information.

Figure 3 shows the spectrum of an observed Type II SN compared to a spectral synthesis model. The comparison demonstrates that current models can be quite successful in reproducing the main properties of observed spectra, and can thus be used to attempt detailed inferrances about composition and structure.Fig. 3. Observed spectrum of SN 2008bk (red), and a model spectrum (blue). Image reproduced with permission from Jerkstrand et al. (2018), copyright by the author(s)

The motivation and goals for spectral synthesis modelling of supernovae and kilonovae are rich and multifaceted. Three of the main science drivers are discussed below.

1) Identifying which elements produce which line features, and determining the elemental abundances in the ejecta. For many years, the challenges to directly infer elemental abundances from supernova spectra appeared almost unsurmountable. Supernovae provide formidable cosmic laboratories with a multitude of complex physics; rapid expansion that Doppler-blends lines, energy cascades over six orders of magnitudes, asymmetries, non-thermal and NLTE effects, molecule and dust formation. Traditional methods from nebular astrophysics run into difficulties and become difficult to meaningfully apply (see McCray 1996, for a good discussion). For the analogous Schuster–Schwarzschild modelling (inner boundary with overlying source-free atmosphere) Mazzali and Lucy (1993)

This led to the curious situation of supernova explosion and nucleosynthesis models having been tested mainly through comparisons to the solar composition, stellar atmospheres abundances, and galactic chemical evolution modelling (e.g., Chiappini et al. 1997; Rauscher et al. 2002; Tominaga et al. 2007), which depend on a long and complex injection series and galactic mixing history of ejecta from a large number of SNe.

Over the last decade or so, however, supernova spectral synthesis models have reached, arguably, a degree of physical realism that this is no longer the case. Synthetic spectra of explosion models are now in quite good overall agreement with observations for several SN classes (see Jerkstrand et al. 2017; Sim 2017, for reviews), and it is possible to attempt to infer abundances to the accuracy needed to directly test individual explosion models and nucleosynthesis theory. By looking into the modelling machinery, we can now also better understand how lines are formed, and from this how to both understand uncertainties and to devise good analytic methods (e.g. Jerkstrand et al. 2012, 2015a, b; Maguire et al. 2018). Thus, one may expect that a paradigm shift is about to occur where this can be systematically done without being limited to more indirect comparisons to estimated stellar atmospheres abundances. There are today clear diagnostic methods established for H, He, C, N, O, Ne, Na, Mg, Si, S, Ar, K, Ca, Ti, Fe, Co, and Ni in supernovae, and direct full-ejecta spectral explosion modelling tests have been done for all major SN classes (e.g. Höflich et al. 1998, 2002; Maeda et al. 2002, 2006; Kasen and Plewa 2005, 2007; Sim et al. 2010, 2012, 2013; Kromer et al. 2010, 2013; Blondin et al. 2013, 2023; Dessart et al. 2013, 2016; Jerkstrand et al. 2012, 2015a, 2016, 2018; Shen et al. 2021; van Baal et al. 2023).

2) Understanding the origin of the elements across the periodic table. With the advent of AT2017gfo, the first kilonova, the upper two thirds of the period table have now also opened up for direct source analysis. While attempts for line identifications and abundance estimates have just begun, diagnostic potential has already been established for _34_Se, _37_Rb, _38_Sr, _39_Y, _52_Te, _57_La, _58_Ce, _60_Nd, and _74_W (Watson et al. 2019; Domoto et al. 2022; Sneppen and Watson 2023; Hotokezaka et al. 2022; Gillanders et al. 2022; Hotokezaka et al. 2023; Pognan et al. 2023). Kilonovae have both advantages and disadvantages when it comes to composition analysis compared to supernovae. The higher velocities ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 0.2$$\end{document} c compared to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 0.02$$\end{document} c) lead to more severe line blending, which complicates identifications. The atomic data for r-process elements is so far much less well known than for elements up to the iron-group, with wavelength uncertainties complicating line matching and A-value and collision strength uncertainties complicating abundance determinations. On the other hand, for KNe there are many more 3D hydrodynamic explosion models with realistic nucleosynthesis available to compute spectra for. The fact that all of the ejecta are radioactive in kilonovae also circumvents a long-standing difficulty of satisfactorily treating the mixing between radioactive regions (^56^Ni-rich) and the rest of the ejecta in supernova modelling.

By systematic spectral modelling of SNe and KNe - event by event, class by class - a major objective is to synthesize a theory for the origin of elements based on direct nycleosynthesis inferrences.

3) Determining the progenitor stellar systems and the explosion mechanisms of supernovae and kilonovae. For collapsing massive stars, the inner parts of the star form a compact object (neutron star or a black hole), and the outer parts are violently ejected. That exotic region of a collapsing stellar core provides a cosmic laboratory where fundamental physics can be probed in regimes not accessible anywhere else in cosmos. The equation of state at the highest densities, strong-field general relativity, neutrino physics, including the enigmatic neutrino oscillations, accretion processes, magneto-hydrodynamics in the context of stellar rotation, and explosive nucleosynthesis are manifested over a few seconds to together produce the supernova phenomenon. Similarly, the explosion of white dwarfs probes a unique physical situation, either arising as two white dwarfs merge, or a single white dwarf accretes matter and initiates collapse as it exceeds the Chandrasekhar limit. By analysing supernova spectra, we can diagnose both the life of the progenitor star, through its hydrostatic nucleosynthesis yields, and its death, through its explosive nucleosynthesis yields and resulting morphology of the ejecta.

In kilonovae, somewhat different high-energy physics is probed. There are here more distinct ejecta components (see e.g. Shibata and Hotokezaka 2019, for a review), ranging from the “contact-squeezed” (probably relatively high electron fraction \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_e=n_p/(n_p+n_n)$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_p$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_n$$\end{document} are proton and neutron abundances, probes the light r-process), “tidal tail” (low \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y_e$$\end{document} , probes heavy r-process), “disc shocks” (probes hypermassive neutron star physics) and “disc viscous” components (probes accretion physics). The study of these outflows opens up a probe of material having been as close to a black hole as matter can get, and can give us insights into the physics of accretion disks and their magnetohydrodynamic processes. The bulk of the ejecta are currently believed to originate from the accretion disk outflows (the last two channels mentioned above). The angular momentum-driven “disc viscous” outflow is particularly important, giving slower ejecta compared to the other channels, as the outflow occurs from the outer edge of the disk where the escape velocity is relatively lower. The lower velocities lead to formation of more narrow emission lines which blend less and therefore give more favorable conditions for identification and analysis. Kilonova spectra give us a unique window on the r-process nucleosynthesis as it happens in real time, and a powerful diagnostic for compact object merger physics.

Sketch of physical situation

Table 1 lists some of the basic physical parameters of CCSNe, TNSNe and KNe. For the phases we will mainly concern ourselves with here, some basic properties of the SN and KN nebulae are

• Homologous expansion. Following an explosion, the faster fragments will soon be further away than the slower ones, and all of them will eventually travel on close to radial trajectories away from the explosion centre. When \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{r} \approx \varvec{v}/t$$\end{document} , we say that the nebula is in homologous expansion. Every point then sees all other points receding away from it, like the galaxies in the Hubble flow. The characteristic velocity scale of the material is around 5000 km \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document} for CCSNe and around for CCSNe, around 8000 km \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document} for TNSNe, and around 50,000 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document} for KNe. KNe have about 100 times less ejecta mass than SNe ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 0.05~M_\odot ~\text{ vs } \sim 5~M_\odot $$\end{document} ), but about the same kinetic energy ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 10^{51}$$\end{document} erg), leading to a factor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 10 higher velocities ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v \propto \sqrt{E/M}$$\end{document} ). Homology is reached faster for more compact progenitors, taking a few hours or days for CCSNe (exploding system of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10-10^3\ R_\odot $$\end{document} ), but only minutes for KNe ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 10^{-4}\ R_\odot $$\end{document} ).
• Radioactive power source. The main power source for the electromagnetic display of many CSSNe and all TNSNe is believed to be the decay of ^56^Ni/^56^Co. Most of the decay occurs in the form of gamma rays, which Compton scatter in the nebula, depositing their energy. Kilonovae are also powered by radioactivity, but here a large number of (r-process) nuclides contribute. The main mediators here are not gamma rays but leptons, and in some cases also \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} -particles and fission fragments. Thus, the power source has here a different distribution with respect to the nebula as a whole, and a different time evolution and thermalization physics. Nevertheless—modelling of both types of nebulae starts with modelling of how radioactivity powers a homologously expanding gas.
• Semi-transparency, with localized line interactions. Supernovae are often discussed and modelled, at different phases, in the theoretical limits of being optically thick (“photospheric”) and optically thin (“nebular”). The true situation is often somewhere in between, and very wavelength-dependent (Jerkstrand 2017). While the optical range clears of continuum opacity relatively quickly, line opacity can give lingering radiative transfer effects for years or decades (Jerkstrand et al. 2011). It has also been understood quite recently that the non-thermal nature of the powering makes UV opacity important, as much energy is reprocessed to UV emission also when temperatures are low (Fransson and Jerkstrand 2015). At all times, much of the optical and IR spectra of both SNe and KNe are formed by fluorescence (Shingles et al. 2023; Pognan et al. 2023).The current picture is well described as that some lines and wavelength ranges can be understood and modelled with quite simple formation physics, e.g. pure scattering (P-Cygni lines) or optically thin emission, whereas others are more complex in formation. The situation is nuanced and much recent work has gone into trying to understand which regimes apply to which lines, and when.The large velocity scale covered by the ejecta (thousands of km/s) implies that line transfer (occurring over the thermal line width, of order 1–10 km/s) occurs over a region small compared to the overall size of the nebula. After a certain spatial distance the comoving (Lagrangian frame) rest wavelength of the photon will have redshifted out of resonance with the transition, and the “Hubble flow” dynamics means it can never come into resonance with that line again, anywhere else. Line interactions are therefore strictly local, and furthermore millions or sometimes billions of resonance scatterings in this local region can be simplified to a a single scattering description in the so called Sobolev formalism (Sobolev 1957; Castor 1970). This simplification has large consequences for convergence properties and algorithm choice, which we will get back to later. Table 1. Characteristic properties of CCSNe, TNSNe, and KNePropertyCCSNTNSNKN \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$M_{ej} (M_\odot )$$\end{document} 2–1510.05 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{ej} (\text{ km } \text{ s}^{-1})$$\end{document} 2000–8000800030,000Composition (Z)1–3020–3030–100 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _{\rm{peak}} (\text{ g } \text{ cm}^{-3})$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-13}$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-12}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-14}$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-13}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-14}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{\rm{peak}} (K)$$\end{document} 5000–20,00010,00010,000 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{\rm{trans}} (d)$$\end{document} 50–20050–10010

Code classes

Codes used for predicting explosive transient light curves and spectra can initially be divided into two groups depending on whether hydrodynamics is solved for or not. The codes including hydrodynamics are referred to as radiation hydrodynamics codes. They often produce bolometric and photometric light curves, but not spectra. Examples include KEPLER (Weaver et al. 1978), the code of Hoeflich and Khokhlov (1996), STELLA (Blinnikov et al. 1998), CRAB (Utrobin 2004), the Bersten code (Bersten et al. 2011), and SNEC (Morozova et al. 2015).

Codes not solving for hydrodynamics must assume a dynamic structure of the ejecta. This structure is, typically, homologous expansion. With the reduction in coding efforts, computing time, and convergence issues related to coupling in hydrodynamics, these codes can instead solve for the radiative transfer to higher accuracy, giving high-resolution spectra. These codes are referred to as radiative transfer codes, and examples include EDDINGTON (Eastman and Pinto 1993), the Höflich code (Hoeflich et al. 1993), PHOENIX (Hauschildt and Baron 1999), CMFGEN (Hillier and Dessart 2012), SEDONA (Kasen et al. 2006), ARTIS/TARDIS (Kromer and Sim 2009; Kerzendorf and Sim 2014), SUMO (Jerkstrand et al. 2011), TH13 (Tanaka and Hotokezaka 2013), SuperNu (Wollaeger et al. 2013), URILIGHT (Wygoda et al. 2019) and POSSIS (Bulla 2023). Many of them do radiative transfer with Monte Carlo methods, especially useful for multidimensional modelling, but for 1D modelling also transfer equation solving can be done. A special subgroup of radiative transfer codes are those considering only local self-absorption in the transfer, through the Sobolev formalism, ignoring the non-local transfer. This ansatz becomes reasonable only at late, nebular times, so these are all NLTE codes. Examples include the 1D codes of Mazzali et al. (2001) and Kozma and Fransson (1998) and the multi-D codes of Maeda et al. (2006), Botyánszki and Kasen (2017),1 and van Baal et al. (2023).

The reason radiation hydrodynamics codes do not do detailed radiative transfer, and vice versa spectral synthesis codes do not do hydrodynamics, has not only to do with issues of coding complexity and convergence; it is the evolving physical situation that sets the needs. Pressure forces, driving dynamics, are important when the density, temperature, and radiation field intensity are high. But those conditions also bring about LTE, broadly speaking, and the internal radiation field which governs the dynamics is well solved for by a simple diffusion approximation. When densities are lower, such that better transfer than diffusive is needed, the ejecta are typically also in their coasting stage (having reached close to the final velocities) and solving for hydrodynamics becomes unnecessary.

The natural division into certain discrete regimes for the physics of explosive transients has both allowed for progress, by having specific efforts tackle specific regimes, but also led to a literature of mostly piece-wise modelling efforts. Ideally, different modellers, specialising in different parts of the problem, would work together and pipeline their efforts. The recent large radiative transfer code comparison project (Blondin et al. 2022) is an important step in this direction, comparing methods and outputs between codes, and identifying paths forward for the field. It continues in the spirit of Blinnikov et al. (1998), who presented the first such efforts for supernova codes.

Review goals and structure

The purpose of this review is to give an overview of the basic physical ingredients for modelling spectra of supernovae and kilonovae, and the numeric implementation of these in modern spectral synthesis codes. Previous related reviews include those of SN spectral formation in the photospheric phases by Sim (2017) and nebular phases by Jerkstrand (2017). This topic is, of course, a vast one and the review only aims to cover the most central themes and aspects. In fact, one of the four central blocks of such modelling—the radiative transfer—is not treated in depth as a section in itself, but instead woven into the other three blocks (temperature, level populations, and radioactive powering). This is partly motivated by that that radiative transfer would truly need its own full review, but also that there exists already many extensive texts on the theoretical foundation and various radiative transfer solution methods (e.g. Baron et al. 1996; Pinto and Eastman 2000; Lucy 2005; Hillier and Dessart 2012; Noebauer and Sim 2019).

We therefore instead go at some depth into the other three blocks, which have received less attention in previous reviews. In Sect. 2 we study the law of energy conservation and its implementation, and how it serves as the most central equation governing the temperature evolution. In Sect. 3 we study how rate equations govern ion and level populations in NLTE modelling, which is now the frontier for both SN and KN modelling. In Sect. 4 we study the radioactivity that powers SNe and KNe, and how its treatment varies between the LTE limit and various types of NLTE.

The target audience of the review are both beginning researchers who are entering the field of computational modelling of astrophysical transients, but also senior SN and KN researchers who have experience with transient analysis and would like to develop a deeper understanding of how models are constructed and can help in interpretation of data. I hope also that the review is useful for the community of experienced modellers for an overview of where the field currently stands. For all groups, my hope is that a mixture of historical context, ingoing physics, numeric formulation, solution techniques, and example illustrations will give the right mixture of aspects to make the review a useful resource.

The following topics are not covered by this review in any significant depth:

• Light curve modelling.
• Comparisons between models and observations.
• Analytic and semi-analytic models.
• Other powering situations than radioactivity.
• Atomic data.
• Numerics of radiative transfer modelling of other explosive transients such as novae, Active Galactic Nuclei, or Tidal Disruption Events.

Temperature

The most fundamental quantity for SN and KN physical modelling is, arguably, the local gas temperature. Temperature governs the ionization and excitation, sets the regime for the spectral formation, and strongly influences the opacity, which in turn governs the radiative transfer. For a model to achieve good fidelity it needs an accurate machinery to compute the gas temperature. It is therefore suitable that we start our survey of the physical modelling of SNe and KNe by looking at how temperature is physically determined, and which equations, approximations, and numeric implementations are in use in current modelling.

Due to the very large cross section for elastic particle collisions at low energies, an equilibrium Maxwell–Boltzmann distribution is typically a good approximation for the bulk of the particles (electrons, atoms and ions), which are at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim $$\end{document} eV energies in SNe/KNe. The quantity 1/2kT is the average kinetic energy per particle, per degree of freedom (three in total from three different spatial directions), in this distribution. This is the fundamental meaning of temperature in this context, and the equilibration is why it is a meaningful and central quantity.

The energy equation

The fundamental physical law governing temperature is energy conservation. The energy conservation equation, in the comoving frame, and ignoring conduction and dissipation which are unimportant in SNe and KNe, can be expressed as (e.g. Hubeny and Mihalas 2014, their Eq. (16.33)):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{\frac{D e}{D t} + p\frac{D}{Dt}\left( \frac{1}{\rho }\right) = \frac{4\pi }{\rho }\int _0^\infty \left( \alpha_\nu J_\nu - \eta _\nu \right) d\nu + \epsilon }}, \end{aligned}$$\end{document}

where D/Dt is the Lagrangian derivative, e is the internal energy per unit mass, p is the gas pressure, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} is the density, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha_\nu $$\end{document} is the (total) absorption coefficient at frequency \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _\nu $$\end{document} is the (total) emission coefficient, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_\nu $$\end{document} is the mean intensity of the radiation field, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon $$\end{document} is a source term (erg \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {g}^{-1}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document} ) representing e.g. energy injection by radioactive decay.

The internal energy e consists of translational kinetic energy by random motions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_k$$\end{document} plus the potential energy of excited/ionized states \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_p$$\end{document} 2 (e.g. Mihalas and Mihalas 1984; Hillier and Dessart 2012):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{ e = e_k + e_p = \frac{3}{2}\frac{kT\left( 1+x_e\right) }{\bar{A}m_p} + \frac{1}{\rho }\sum _{\rm{ion,exc}} n_i E_i}} \,, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_e \equiv n_e/n_{nuclei}$$\end{document} is the electron fraction, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{A}$$\end{document} is the mean atomic weight (number), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_i$$\end{document} is the number density of state i, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_i$$\end{document} is the energy of that state. For atoms and ions “excitation” involves change of electronic orbitals, whereas for molecules the nuclei can also be excited into rotational and vibrational states.

Absorption and emission can occur by a variety of processes, which can be divided into the categories of free-free, free-bound/bound-free, and bound-bound, referring to the initial and final state of the electron involved in the process. Another division axis is along scattering and absorption, referring to whether photons merely change direction in the interaction process3 or are destroyed. Because scattering does not change the states of the material particles, if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi _\nu $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _\nu $$\end{document} have explicit scattering components, these may be removed before use in Eq. (1). Temperature generally has a strong influence on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\chi _\nu $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta _\nu $$\end{document} which, together with the explicit T-dependency of e (Eq. (2)), is why we can say that the energy equation is the key equation for setting the temperature. For example, the contribution to the emission integral by a line ul in LTE from an ion i is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4\pi /\rho \int j_\nu d\nu = (n_i/\rho ) A_{ul} g_u \exp {\left( -E_u/kT\right) }/Z(T)$$\end{document} , where Z(T) is the partition function.

When NLTE is considered4, one may physically equivalently, but computationally sometimes more expedient and numerically better conditioned, write an equation describing the evolution of the thermal kinetic energy stored in the pool of atoms, ions and thermal electrons (e.g. Kozma and Fransson 1998; Jerkstrand 2011):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{\frac{De_k}{Dt} + p\frac{D}{Dt}\left( \frac{1}{\rho }\right) = h - c \,,}} \end{aligned}$$\end{document}

where h is the heating per unit mass (energy flow going towards increasing the kinetic energy of particles) and c the radiative cooling per unit mass (energy flow going towards decreasing the kinetic energy of the particles). Adiabatic cooling is represented by the second term on the LHS. For example, if 50% of the radioactive decay power (r.p) ends up as internal kinetic energy of the particles, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{r.p.}=0.5 \epsilon $$\end{document} . As another example, the cooling through a line transition is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c_{lu} = \rho^{-1} E_{ul}\left( q_{lu}(T)n_e n_l - q_{ul}(T) n_e n_u E\right) $$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{lu}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{ul}$$\end{document} are upward and downward collision rates, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_e$$\end{document} is the electron number density.

Note that charge exchange reactions and chemical reactions (molecule and dust formation and destruction) are also associated with heating and cooling which needs to be accounted for if these processes are included. The reader is referred to Hillier (2003) for further discussion of the choice of energy equation.

The thermal equilibrium approximation corresponds to ignoring the two time-derivative terms in the energy equation (Eq. (1) or, equivalently, Eq. (3)), so the whole LHS becomes zero. This changes the equation from a first-order initial value problem into an algebraic equation. If any source term entries ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon $$\end{document} ) are considered as “radiation” (e.g. the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} particles from radioactive decay), or if the source term is zero, thermal equilibrium is equivalently sometimes called radiative equilibrium, as (from Eq. (1)) the absorption of radiation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\epsilon + 4\pi /\mathbf {\rho } \int \alpha_\nu J_\nu d\nu )$$\end{document} balances the emission of radiation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(4\pi /\mathbf {\rho } \int \eta _\nu d\nu )$$\end{document} .

For the case of Eq. (3), the thermal equilibrium condition is simply

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{h = c.}} \end{aligned}$$\end{document}

The expression “algebraic equation” is here to be taken as a heuristic description. Both heating and cooling terms are, in practice, implicitly obtained by solving associated large numeric problems. For example, heating may be determined by a full Monte Carlo simulation of the gamma-ray Compton scattering process, and cooling by solving thousands of linked rate equations.

When is the thermal equilibrium approximation suitable? It holds when the time-scales of both heating and cooling are fast compared to the time-scales over which the fundamental physical situation—density and powering—changes. The latter are characterised by the expansion and power source time-scales, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho /\dot{\rho }$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon /\dot{\epsilon }$$\end{document} , respectively. For the case of homologous expansion and single-isotope radioactive power these expressions equal \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho /\dot{\rho }=t/3$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon /\dot{\epsilon }=\tau _{\rm{decay}}$$\end{document} . Section 4.2 in Jerkstrand (2011) contains a further discussion and some illustration of the accuracy of the approximation for SNe, and Pognan et al. (2022b) makes a study of it for KN applications. For the KN case one can show that, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^{-1.3}$$\end{document} r-process raw radioactive decay power and a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( 1+t/t_b\right) ^{-1.5}$$\end{document} thermalization factor (Kasen and Barnes 2019; Waxman et al. 2019, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_b$$\end{document} is a thermalization timescale), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon /\dot{\epsilon }$$\end{document} becomes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t/\left( 1.3 + 1.5\left( 1+t/t_b\right) ^{-1}\times t/t_b\right) $$\end{document} . This expression limits to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t/1.3\ (=0.77t)$$\end{document} for small t and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t/2.8\ (=0.36t)$$\end{document} for large t. Thus, the homologous time-scale t/3 always puts more stringent constraints for KNe.

Several radiative transfer codes, e.g. STELLA, SuperNu, CMFGEN, and (since 2022) SUMO have the capacity to retain the time-derivative terms in the energy equation (STELLA and SuperNu use Eq. (1), SUMO uses Eq. (3), whereas CMFGEN checks that both Eq. (1) and Eq. (3) are fulfilled)—and thus can avoid the issue of the accuracy of the thermal equilibrium approximation. Numerically, implicit time-differencing gives equations of the same form as when doing thermal equilibrium solutions; the implicit time derivative just adds in more source terms, known from the previous time step. Thus, doing time-dependent solutions has from the implementation point-of-view no real additional complexity or convergence issues compared to equilibrium modelling. The difference is only the need to evolve the problem as an IVP through many time steps. If just a single epoch is of interest, having to do a large number of previous epochs is then the “cost” of dropping the approximation. But if one desires a full time series of spectra anyway, there is no real extra cost.

Solving the time-dependent energy equation is done by a finite difference scheme. STELLA and CMFGEN use a fully implicit time discretization (all terms defining the derivative evaluated at the new time step, Blinnikov et al. 1998; Hillier and Dessart 2012), whereas SuperNu applies a semi-implicit one (some terms defined at the new time step, some at the previous, Wollaeger et al. 2013). The performance of fully implicit versus semi-implicit schemes can often vary significantly with application (e.g, Stone et al. 1992). When applied to a single ODE with constant coefficients, textbooks show us that a fully implicit differencing scheme is unconditionally stable and first order accurate (meaning that accuracy is linearly proportional to the numeric step size), whereas a semi-implicit scheme is unconditionally stable but second order accurate (meaning that accuracy is quadratically dependent on numeric step size, which is the highest order achievable for unconditional stability). However, these statements do not hold universally for all problems. Here we are dealing with an ODE with non-constant coefficients, and in addition this equation is typically co-solved with other non-linear equations. As discussed by Press et al. (1992), implicit methods are then usually, but not always, stable, and accuracy cannot be so precisely pre-stated in general. Hoeflich et al. (1993) parametrizes the degree of implicitness, but states that in practice the fully implicit limit is used for all model runs. Lucy (2005) analyses accuracy improvements one may obtain by using information from two previous timesteps instead of one, for the case of solving the energy equation for the radiation field energy, and the first-order moment equation for radiative transfer.

When a time-dependent energy equation is used, there are normally also other equations in the system with time-derivative terms retained (giving a set of coupled ODEs). The choice of time-step then comes down to which of these equations requires the shortest step for stability and accuracy. In addition, the task of identifying an optimal discretization method (among implicit, semi-implicit, Runge–Kutta, Richardson extrapolation, predictor-corrector, etc.) becomes very complex. In the SN/KN literature so far, the most common approach is therefore to keep it simple with an implicit method using a fixed, moderate-sized time-step, typically around 10%. This time step ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t=0.1t$$\end{document} ) is always shorter than the homologous expansion time-scale (0.33t) meaning that time resolution with respect to density changes is controlled. The other time-scale to keep track of is the power source time-scale, which is the decay time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{\rm{decay}}$$\end{document} in the case of single-isotope radioactivity (111 d for ^56^Co). Desiring to stay in the regime \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t \le $$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{\rm{decay}}/\left( \text{ factor } \text{ few }\right) $$\end{document} then means that 10% stepping is sufficient for epochs up to a few years. Conveniently, this is also when other more long-lived isotopes like ^57^Co ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{\rm{decay}}=392$$\end{document}  d) and ^44^Ti ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{\rm{decay}}=87$$\end{document}  y) take over in SNe. For KNe, we can estimate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau _{\rm{decay}} = (0.36-0.77)t$$\end{document} as discussed above. Thus, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 10% time-stepping has a quite solid physical anchoring for both SN and KN applications, for any conceivable range of epochs.

A final note is, however, that having \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta t \ll \left[ 0.33t,\tau _{\rm{decay}}\right] $$\end{document} is a first, rough requisite; what really matters is to well resolve the evolution of the physical state variables. Should these vary with time more rapidly (e.g., a sudden episode of dust formation that significantly changes the cooling function), a yet better time-resolution may be needed. Ideally is an adaptive time-step control implemented which can monitor this and redo the time-step if a violation occurs.

Heating

When working with e (energy Eq. 1), the term “heating” refers to processes that increase the internal energy, whether it is kinetic energy of the material particles, or excitation/ionization energy. When working with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_k$$\end{document} (energy Eq. 3), it refers to processes in the first category only.

In layers exposed to radioactive decay particles (see Sect. 4 for details), heating by impact and cascading of these will be important. Heating by absorption of the diffuse radiation (photoelectric absorption, free-free absorption, and bound-bound absorption followed by collisional deexcitation) can, however, also play a role. Thus, heating should be computed as the sum of a radioactive heating estimate and a radiation field heating estimate.

NLTE

Radioactive heating. In NLTE one needs to compute the fractions of radioactive decay energy that goes to random thermal motions, and which go to excitations and ionizations. However, for the purposes of the energy equation, as long as the gas is at least singly ionized or so (Kozma and Fransson 1992) no big error is made if one takes the full deposition power to go to heating. In SNe, the ionization decreases with time though, so this approximation gets progressively worse and beyond some point in time the fraction needs to be computed—more details on this in Sect. 4. For KNe, it has recently been understood that ionization turns around and increases again after the diffusion phase ends (Hotokezaka et al. 2021; Pognan et al. 2022b). The ionization state is always high enough that no big error is incurred if the 100% heating approximation is retained at all times for the energy equation.

Radiation field heating. By whichever radiative transfer method is used, radiation field heating by bound-free and free-free processes are straightforwardly computed. For bound-bound, line absorption generate photoexcitation which lead to increased upper level population in the next iteration. This in turn modifies the net bound-bound cooling rate for the electrons on transitions involving that level (Eq. (3)). In an average sense, the level push-ups induced by photoexcitations will lead to a smaller net cooling as more downward (heating) paths become available and fewer upward (cooling) ones. Thus, here line absorptions do not enter explicitly as a heating term in the energy equation, but instead they will degrade the cooling term. This is what happens in NLTE codes like SUMO and CMFGEN.

LTE

Radioactive heating. In LTE modelling one works with Eq. (1), and the radioactive heating is just the total locally deposited energy, with no need to specify its components.

Radiation field heating. For the radiation field heating estimate, Monte Carlo codes typically use the volume accumulator method of Lucy (1999) to estimate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_\nu $$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} J_\nu = \frac{1}{4\pi \Delta \nu V}\sum _{\Delta \nu } \epsilon \times \Delta l, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \nu $$\end{document} is the frequency bin width, V is the cell volume, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon $$\end{document} is the packet energy, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta l$$\end{document} the travel segment. To compute a radiative heating rate from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_\nu $$\end{document} , LTE codes require specification of what fraction of absorbed radiation thermalizes and what fraction scatters, i.e. how the total absorption coefficient breaks down into absorption (thermalization) ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _{\mathrm{abs,\nu }}$$\end{document} ) and scattering ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _{\mathrm{scatt,\nu }}$$\end{document} ). The radiative heating depends only on the absorption part:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{h_{\mathrm{rad-field}}^{\rm{LTE}} = \rho^{-1} \int _0^\infty 4\pi J_\nu \alpha _{\mathrm{abs,\nu }} d\nu }}. \end{aligned}$$\end{document}

While this division is straightforward for continuum processes, SNe and KNe tend to become dominated by lines for the opacity. Determining what happens following a line absorption is difficult and computationally demanding - and it turns out that the importance of fluorescence clashes with the omission or approximative description of this process in many LTE codes - more on this in Sect. 2.3.

Cooling

When working with e (energy Eq. 1), the term “cooling” refers to processes that decrease the internal energy, whether it is kinetic energy of the material particles or excitation/ionization energy. When working with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_k$$\end{document} (energy Eq. 3), it refers to processes in the first category only.

Cooling will occur by various radiative processes, of which line cooling is typically the most important in SNe and KNe. Other cooling processes include recombination and free-free emission, which however typically contribute by of order a few percent only.

NLTE

The line cooling by electrons (used for Eq. (3)) is given by:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} c = \rho^{-1} h\nu _0 n_e \left( q_{up}(T) n_l - q_{\rm{down}}(T)n_u\right) . \end{aligned}$$\end{document}

The first thing to note is that if one inserts the exact LTE values for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_l$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_u$$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_u/n_l = g_u/g_l \exp (-E_{ul}/kT)$$\end{document} )), c becomes identical zero, because \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{up} = q_{\rm{down}} g_u/g_l \exp {\left( -E_{ul}/kT\right) }$$\end{document} . But the upper level still emits photons at a rate of A (or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A \times \beta _s$$\end{document} , more below), which must represent (at least partially) cooling. In fact, it is this radiative leakage that gives the upper level a population slightly below the LTE limit that makes all the difference. That deviation from LTE is proportional to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1/n_e$$\end{document} in magnitude if collisions dominate. Thus, as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_e$$\end{document} becomes large, Eq. (7) shows that the cooling is the small difference between between two large terms. Clearly care both in formulation and numerics is needed here.

Calculating cooling by the direct summing of level populations times A or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A\times \beta _s$$\end{document} (instead of using Eq. (7)) can be done only in thermal excitation-only NLTE making a consistent choice in the rate equations, as in e.g. Hotokezaka et al. (2021). However, as soon as other processes are allowed to affect the NLTE solutions (e.g. recombination), one instead has to work with the electron net collision rates.

In general, how much a given atom or ion can cool (by being a target of collisional excitation by the thermal electrons) depends on its NLTE state, which in turn depends on both the number density of electrons, the number density of the parent ion (affects recombination inflows), and the radiation field (affects photoexcitation and photoionization rates). However, in the low-density limit most atoms/ions are in the ground state, and one can assume that each collisional excitation leads to radiative decay (cooling). There then exists a unique cooling function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document} that depends on T only, so that cooling by ion i is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} c_i = \rho^{-1}\Lambda _i(T) n_i n_e \,. \end{aligned}$$\end{document}

Figure 4 shows the low-density cooling functions of Nd II, as well as of different ions of Ce. The rapid growth with temperature is a characteristic feature of cooling functions due to the exponential terms that enter expressions for the collisional excitation rates. It is what keeps the temperature range obtained in astrophysical nebulae quite limited, with H II regions, planetary nebulae, SN and KN ejecta all being in the range \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 3000–30,000 K for the most part. The Nd II curve illustrates that both forbidden (red, dotted) and allowed (blue, dashed) transitions can be important to include in the cooling modelling. The Ce curves illustrate that cooling can be sensitive to the ionization state of the gas, so there is an intrinsic strong coupling betweeen temperature and ionization.Fig. 4. Image reproduced with permissions Left: Cooling function of Nd II in the low-density limit, from Hotokezaka et al. (2021). Right: Cooling functions of different ions of Ce, in the low density limit, from Pognan et al. (2022b). Images reproduced with permission, copyright by the author(s)

LTE

In LTE, a subtlety affects the formulation of the line cooling. One plausible approach would to take the line cooling as the Boltzmann level populations times \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A\beta _s$$\end{document} (and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h\nu _0)$$\end{document} , which is the effective decay rate. Another plausible approach is to set the source function equal to the Planck function. One may show that this corresponds to taking the Boltzmann level populations times A. The source if this discrepancy is that the LTE approach at its very basic level lacks self-consistency in its formation (see e.g. discussion in Chapt. 2 of Mihalas 1978). Here, self-trappings change the radiation field from a Planck function, and the level populations from pure LTE, but that cannot be accommodated for in the framework. It is also significantly easier/faster to use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_\nu = B_\nu $$\end{document} , and one may argue that engaging in other more expensive treatments has little meaning as yet further self-inconsistencies are added. However, the Kirchoff–Planck treatment may overestimate the cooling (and thus underestimate temperatures) when optically thick lines are important.

With the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_\nu = B_\nu $$\end{document} approach, cooling in LTE is computed as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{ c^{\rm{LTE}} = \rho^{-1}4\pi \int _0^{\infty } \alpha _{\text{abs},\nu } B(\nu ,T) d\nu }}, \end{aligned}$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{\alpha _{\mathrm{abs,\nu }} =\alpha ^{\text{abs}}_{bb,\nu } + \alpha _{bf,\nu } + \alpha _{ff,\nu }}}. \end{aligned}$$\end{document}

The second and third terms on the RHS are due to continuum heating (bound-free and free-free, respectively), which have well-known expressions. The first term on the RHS represents absorption in lines, meaning photoexcitation followed by collisional de-excitation. The “true” absorption (“thermalization”) probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^{\text{abs,true}}_{bb}$$\end{document} is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{ p^{\rm{abs,true}}_{bb} = \frac{\sum _{\rm{down}} q_{\rm{down}}n_e}{\sum _{\rm{down}} \left[ q_{\rm{down}}n_e + A\beta _s\right] } \,, }} \end{aligned}$$\end{document}

which is typically small, in the range \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-6}-10^{-4}$$\end{document} as long as allowed de-excitation transitions ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A\gtrsim 10^6\ \text{ s}^{-1}$$\end{document} ) exist (Kasen et al. 2006).

However, it is expensive to compute \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{bb}^{\text{abs,true}}$$\end{document} for all lines. Many LTE codes therefore use a parameterized constant value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{bb}^{\text{abs}}=\epsilon _0$$\end{document} , letting the photon thermalize with probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _0$$\end{document} and resonance scatter with probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1-\epsilon _0$$\end{document} . An additional advantage of this is that lines can be binned and added up to an “expansion opacity”, as what follows after a line absorption becomes the same for all of them, and therefore the order between the lines in a bin does not matter. Such an expansion opacity is computed as, for wavelength bin i,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \alpha _i = \frac{1}{ct} \frac{1}{\Delta \lambda _i} \sum _{j} \lambda _j \epsilon _0 \left( 1-\exp {\left( -\tau _{s,j}\right) }\right) . \end{aligned}$$\end{document}

Surprisingly, the value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _0$$\end{document} must be high ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 1 $$\end{document} ) for realistic light curves and spectra to be produced (Blinnikov et al. 1998; Kasen et al. 2006; Kozyreva et al. 2020). This is interpreted to be because the lack of fluorescence in this treatment gives an artificial lack of wavelength transformations for the photons. A high \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _0$$\end{document} value is a way to heuristically mimic the fluorescence process, but this comes at the cost of introducing errors in the energy equation, as more radiative energy is thermalized than it should be.

Improvement on the situation can be done by treating at least some of the lines in detail, including the fluorescence. This is done by SEDONA (Kasen et al. 2006), which splits

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \alpha _{bb,\nu }^{\rm{abs}} = \alpha _{bb,\nu }^{\rm{direct}} + \alpha _{bb,\nu }^{\epsilon _0} = \chi _{bb,\nu }\times \left( f_{\rm{direct}}p_{bb}^{\rm{abs,true}} + (1-f_{\rm{direct}})\epsilon _0\right) . \end{aligned}$$\end{document}

Using a thermalization probability \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{bb}^{\rm{abs,true}} \ll 1$$\end{document} for some lines and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 1$$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _0$$\end{document} ) for others means that the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _0$$\end{document} -treated ones will be most influential for determining the temperature. In LTE with thermal equilibrium, and assuming that lines are dominant for the opacity, we have a condition

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{ h_{\rm{radioactivity}} + \int _0^\infty \frac{4 \pi }{\rho } J_\nu \alpha _{bb,\nu }^{\rm{abs}} d\nu = \int _0^\infty \frac{4\pi }{\rho } B_\nu (T) \alpha _{bb,\nu }^{\rm{abs}} d\nu . }} \end{aligned}$$\end{document}

If heating is radiatively dominated (with the microphysically correct thermalization probability), we see that the scale of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _{bb}^{\text{abs}}$$\end{document} does not actually matter—we can multiply this function with any number we want and the LHS and RHS of Eq. (14) would change by the same factor in the limit that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{\rm{radioactivity}}$$\end{document} is unimportant. LTE codes with expansion opacities therefore obtain temperatures quite robust to this uncertainty in the photospheric phase, especially in the layers outside of the radioactive material, which are also the most important ones for the spectral formation.

When radioactivity dominates the heating, however, we have (for a roughly gray \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _{bb}^{\text{abs}}$$\end{document} ):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {{ aT^4 \approx \frac{h_{\rm{radioactivity}}}{\alpha _{bb}^{\text{abs}}}. }} \end{aligned}$$\end{document}

If \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _{bb}^{\text{abs}}$$\end{document} is too large by a factor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^4$$\end{document} (using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon _0 \sim \ 1$$\end{document} instead of the microphysically accurate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{bb}^{\text{abs,true}} \lesssim 10^{-4}$$\end{document} ), the temperature estimate will therefore be a factor 10 too low. This is likely one driving factor (combined with deviation from LTE for level populations) for why LTE codes give significantly lower SN and KN temperatures than NLTE ones from a quite early phases (e.g. Pognan et al. 2023; Blondin et al. 2023). One can also note that KNe are quite different with respect to SNe in that here all the ejecta layers are radioactive. This raises also the possibility that LTE KN temperatures in the photospheric phases may be less accurate than in SN modelling.

Temperature from radiation field

A gray model for a (source-free) plane-parallel stellar atmosphere gives a temperature stratification going from \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_{\rm{eff}}$$\end{document} at the photosphere to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.81 T_{\rm{eff}}$$\end{document} at the surface (Mihalas 1970). From this, the choice of temperature in an isothermal model is often \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_e=0.9 T_{\rm{eff}}$$\end{document} , an average value for the layer. Non-gray models tend to give somewhat lower atmospheric temperatures so the coefficient can also be somewhat lower (e.g. 0.7 in Lucy and Solomon 1970).

In diffusion-phase SNe and KNe a photosphere can analogously be defined from where the optical depth to the surface is 2/3, although certain issues become more severe here compared to the stellar atmospheres case (stronger frequency-dependency of the photospheric depth, and a bigger differences between thermalization depth and photosphere, see Sim 2017, for a good discussion). For the analogous Schuster–Schwarzschild modelling (inner boundary with overlying source-free atmosphere) Mazzali and Lucy (1993) proposed that for SNe a more suitable approximation is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_e = 0.9 T_R$$\end{document} , where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} T_R \equiv \frac{\langle h \nu \rangle}{3.82k}. \end{aligned}$$\end{document}

Here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\langle h \nu \rangle$$\end{document} is the mean energy of photons in the radiation field, obtained from either a radiative transfer solution or a Monte Carlo simulation. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_R$$\end{document} can be said to be the temperature of a blackbody with the same “color” as the radiation field.

Another variant was proposed by Kromer and Sim (2009), who in their “simple temperature/ionization” methodology use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_e = T_J$$\end{document} , with

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} T_J \equiv \left( \frac{\pi J}{\sigma }\right) ^{1/4}, \end{aligned}$$\end{document}

which is the temperature of a blackbody field with the same intensity as the radiation field. Here J is the bolometric radiation field strength and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} is the Stefan–Boltzmann constant. As the radiation field gets more dilute with time, Eq. (17) will inevitably start to underestimate the temperature even if gas-photon coupling is still strong.

Kromer and Sim (2009) compare results of this method with a more accurate (named “detailed temperature/ionization”) method in which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_e$$\end{document} is instead solved from balancing heating with cooling as described above (Fig. 5). For the radiative heating estimate, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4\pi \int \alpha_\nu J_\nu d\nu $$\end{document} , in this detailed approach, the authors implement the method of Mazzali and Lucy (1993) for a 2-parameter model of the radiation field (there used to estimate ionization rates):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} J_\nu = W B_\nu (T_R) \,, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$W \equiv J/B(T_R)$$\end{document} is called the dilution factor. Thus, here the requirement of Eq. (17) that the radiation field should be as strong as a blackbody is relaxed, while the requirement of a similar SED shape is retained. The strong coupling assumption is also relaxed as with the thermal equilibrum the radiation field “temperature” can be different than the gas temperature.

Due to their simplicity and particular ease of implementation in Monte Carlo codes, Eqs. (16) and (17) have been quite widely used in the modelling literature, not least for extensive parameter space modelling of photospheric SN spectra using the codes of Mazzali and Lucy (1993), Mazzali (2000), Kerzendorf and Sim (2014). Moreover, they have also carried over to be used for many currently available 3D, full-ejecta models of kilonovae (e.g. Tanaka and Hotokezaka 2013; Bulla 2023).

Figure 5 shows the Kromer and Sim (2009) comparison of temperatures obtained through 17, and from thermal equilibrium, for a Type Ia model. The agreement is quite good for 0–15 d, whereas after this Eq. (17) gives an underestimate by a factor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 2, due to its neglect of the dilution as discussed above. We can also see that using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_R$$\end{document} instead of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_J$$\end{document} to estimate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_e$$\end{document} tends to give a better approximation at all phases. While this gives some information about what errors the simplifications of Eqs. (16) and (17) incur, for one particular SN type (Type Ia), it is currently not known what the level of accuracy is for CCSNe or KNe.Fig. 5Left: Comparison between temperature time-evolution at velocity coordinate 9600 km \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document} in Type Ia LTE SN model using different treatments. Electron temperatures \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_e$$\end{document} using Eq. (17) are plotted as green pluses, while those from thermal equilibrium are plotted as black points. The blue and red points are \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_R$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T_J$$\end{document} in the thermal equilibrium model. Right: Same as the left figure, but showing profiles (temperature versus velocity) at 31 days. The black/violet/blue/green/red/orange vertical lines show the mean radii of last scattering in the U, B, V, R, I bands. Images reproduced with permission from Kromer and Sim (2009), copyright by the authors

Numeric aspects and some comparisons

Energy conservation is defined by a single equation (Eq. (1) or equivalently Eq. (3)) but this contains many terms, most of which are highly non-linear in temperature. In NLTE, the constituent terms are not explicit functions, but evaluated at different temperatures by solving the associated rate equations (described in the Sect. 3). Newton–Raphson is a standard algorithm choice, but for single non-linear equations there are also other options (see e.g. Chapt. 9 in Press et al. 1992). Newton–Raphson is in particular to prefer if there is a good starting guess, and when it comes to temperature it is relatively straight-forward to roughly pre-estimate its value. The implicit nature of the problem in NLTE means the derivative is computed numerically (i.e., the heating and cooling terms typically do not have explicit, manageable analytic expressions, and therefore cannot yield analytic temperature derivatives).

Stability is enhanced by both damping and capping steps. In the thermal equilibrium iterations described by Lucy (1999, 2003), a step damping factor of 0.5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-$$\end{document} 0.8 and a step cap of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\Delta T|<200$$\end{document} K per iteration are reported to be beneficial for robust convergence. In addition, an adaptive scheme is stated to be needed that allows a further reduction of both these in (rare) cases where corrections do not evolve healthily. Lucy (2003) uses a temperature convergence criterium of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta T<$$\end{document} 1 K, and also checks that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(h-c)/c$$\end{document} is small.

It is beneficial for convergence to solve the (single) energy equilibrium equation together with the level population equations (e.g. Lucy 2003), or at least some of them if the total number is too high. The thermal equation increases the system size by only one, while allowing for sometimes significantly fewer iterations. The coupling becomes in particular powerful if the main cooling species are identified and co-solved for.

We close this section by showing a comparison of temperature profiles computed with a variety of codes for a Type Ia SN test model (Fig. 6, from Blondin et al. 2023). As can be seen, temperature predictions for even a simple test model (this one contains only Fe, Co, Ni) vary quite significantly between codes. All models shown here compute temperature from the first law of thermodynamics, but clearly estimates of heating and cooling terms vary quite significantly. The comparison illustrates the challenges for temperature modelling even for state-of-the-art radiative transfer tools, and much work remains to be done to arrive at robust, accurate determination of this fundamental quantity in SNe and KNe.Fig. 6. Temperature profiles obtained by various codes for a Type Ia test model at 40 d (left) and 200 d (right). Image reproduced with permission from Blondin et al. (2022), copyright by the author(s)

Rate equations

We now move on to study numeric techniques for solving rate equations, needed for NLTE modelling. LTE codes make use of the Boltzmann equations for excitation which for given values of T and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} (and composition) provide any level population by a formula depending only on the energy and statistical weight of the level and the partition function, well known quantities. LTE ionization is computed from the the Saha–Boltzmann equations, where also the electron density enters. This makes the equations non-linear, but their solution is a straightforward procedure which will not be discussed further here.

After some time populations start to deviate significantly from LTE. Figure 7 shows computed departure coefficients (NLTE level populations relative to LTE ones) for Ni II in a Type Ia model at 330 d (top), and for Ce II in a KN model at 20 d (bottom). We can see how NLTE effects are strong in both cases, and increase with excitation energy. A few take-away points are:

1. The NLTE populations are in general very different from their LTE ones. This tells us that LTE modelling beyond some point in time is not expected to give even approximately correct results; NLTE gives a completely different state of the material.
2. Departure coefficients initially decline with increasing level number, driven by the fact that there are more and more possible radiative de-excitation paths opening. However, a trend reversal can arise due to recombination.
3. The Ce plot (bottom panel) shows that with enough levels included, the levels fall into two distinct groups; those with relatively high level populations and departure coefficients of order \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-3}$$\end{document} (level numbers \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim 50$$\end{document} in this example) and those with much lower populations and departure coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-10}$$\end{document} – \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-6}$$\end{document} (level numbers \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gtrsim 50$$\end{document} in this example). This split arises as levels over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 50 in Ce III start to have allowed-transition de-excitation paths, whereas the lower lying ones have only forbidden-transition ones. This situation is not to be confused with that the lower group is necessarily more important for emission; it is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n\times A$$\end{document} that governs emission, and the upper group is depopulated because of their large A values. The ionization structure will in general also deviate strongly from LTE, in addition to the excitation structure devitations shown here. This is because ionization occurs mainly by non-thermal electrons, and recombination by radiative and di-electronic processes, both different to the thermal processes assumed to dominate in LTE.Fig. 7Top: Departure coefficients for the first 17 levels of Ni II in a Type Ia model at 330 d (the two types of symbols referes to using two different sources of atomic data). Bottom: LTE level population solutions for Ce II, compared to NLTE values, in a KN model at 20d. “Limited NLTE” means NLTE excluding photoexcitation and photodeexcitation rates. The upper panel shows the number densities, and the lower panel the departure coefficients. Images reproduced with permission from [top] Shingles et al. (2020) and [bottom] Pognan et al. (2022a), copyright by the author(s)

Formulation

We will focus on NLTE rate equation systems for the internal states of atoms and ions in steady state, which means that balance holds between inflow and outflow. For non-steady state effects, the reader is referred to discussions in Fransson and Kozma (1993, late evolution of SN 1987A), Utrobin and Chugai (2005), Dessart and Hillier (2008, photospheric phase of Type II SNe), Fransson and Jerkstrand (2015, nebular phase of Type Ia SNe), and Pognan et al. 2022b, nebular phase of kilonovae). Discussion and examples are taken for atomic modelling; works on solving rate equations for molecules can be found in Liu et al. (1992), Liu and Dalgarno (1995), Liljegren et al. (2020, 2023), Rho et al. (2021).

There are many rate equations for each grid cell, one for each included level. As there may be dozens of elements, maybe half-a-dozen ionization stages, and hundreds or even thousands or levels per ion, the count quickly gets large, often to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \gtrsim 10^5$$\end{document} per cell. Normally one splits the full equation system up into blocks, for example may each ion (e.g. Fe II) define a separate block which then has \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=10^2-10^3$$\end{document} equations. These blocks are then solved sequentially and iteratively. This split and iteration is done because matrix storage scales with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^2$$\end{document} and inversion time with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^3$$\end{document} . Thus, 100 equations would cost \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto 10^4$$\end{document} in storage if coupled, but \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto 10 \times 10^2 = 10^3$$\end{document} if split into blocks of ten, similarly the number of operations would be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto 10^6$$\end{document} if coupled but \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto N_{iter} 10^3$$\end{document} if split into ten blocks; clearly if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_{iter}$$\end{document} (number of iterations needed) is not too large a saving is made. When solving for each ion, all other quantities (populations of levels in other elements and ions, temperature, electron density) are held fixed.

Writing down a balance between inflow (LHS) and outflow (RHS) for level i in a N-level block gives an equation of form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sum _{j=1, \ne i}^N n_j R_{j,i} + s_{i,\mathrm {ext}} = n_i \times \left[ \sum _{j=1, \ne i}^N R_{i,j} + R_{i,\mathrm {ext}}\right] , \end{aligned}$$\end{document}

Here \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{i,j}$$\end{document} are reaction rates (unit \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document} ) between levels i and j, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{i,\mathrm {ext}}$$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^{-3}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document} ) is a block-external inflow rate (a source term) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{i,\mathrm {\text{ext}}}$$\end{document} is a reaction rate ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document} ) to block-external states (a sink term). For single-ion blocks, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{i,\mathrm {ext}}$$\end{document} could for example contain contributions by recombination inflows, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{i,\mathrm {ext}}$$\end{document} could contain contributions by photoionization outflows.

The components of the different terms arise from a variety of collisional and radiative processes, see e.g. Chapt. 4 in Jerkstrand (2011), for an overview. In the optically thin limit, no \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{i,j}$$\end{document} terms depend on any of the level populations in the block so the system is linear, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{R}\textbf{n} = \mathbf {s_{\rm{ext}}}$$\end{document} , where the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{R}$$\end{document} matrix contains the R terms and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbf {s_{\rm{ext}}}$$\end{document} array the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$s_{i,\mathrm {ext}}$$\end{document} terms. If line optical depths are included and treated in the Sobolev formalism (Sobolev 1957; Castor 1970), the effective spontaneous radiative decay rate from level i to j is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$A\times \beta ^s_{i,j}$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta ^s_{i,j}$$\end{document} is the Sobolev escape probability, which for homologous expansion is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \beta ^s_{i,j} = \frac{1- \exp {\left( -\tau ^s_{i,j}\right) }}{\tau ^s_{i,j}} \,, \end{aligned}$$\end{document}

where the Sobolev optical depth \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau ^s_{i,j}$$\end{document} is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \tau ^s_{i,j} = \frac{1}{8\pi }\frac{g_i}{g_j} A_{i,j} \lambda _{i,j}^3 n_j t \left( 1 - \frac{g_j}{g_i}\frac{n_i}{n_j}\right) , \end{aligned}$$\end{document}

which introduces a non-linearity as the R component for the spontaneous decay ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{i,j}^{\rm{spont}}=A_{i,j}\times \beta _{i,j}^s$$\end{document} ) now depends on the level populations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_j$$\end{document} themselves. Another source of non-linearity can arise if one considers also continuum destruction probabilities. If photons are allowed to be destroyed by photoionization on e.g. a level k in the block before escaping the Sobolev resonance, another term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _c(n_k)$$\end{document} would enter (Hummer and Rybicki 1985; Chugai 1987; Jerkstrand 2011).

Solution with Newton–Raphson method

The non-linearity introduced by self-absorption discussed above is dealt with quite effectively by fixed-point iteration (compute \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} values for given populations, hold fixed, solve linear system, repeat, Lucy 1991), but one may also choose to solve for the next corrective step in a multi-dimensional Newton–Raphson scheme,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{J} \cdot \Delta \textbf{n} = -\textbf{e} \,, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{J}$$\end{document} is the Jacobian of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{R}$$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\times N$$\end{document} array of all partial derivatives) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{e}$$\end{document} is the length-N array of current errors of the rate equations (residuals of total inflow minus total outflow for each level, Eq. (19)). To minimize numeric errors, and save computation time, it is beneficial to derive analytic Jacobian entries when possible (while this is typically not doable for the energy equation (Sect. 2), it is manageble for each rate equation, but on the other hand there is a very large number of them). If doing numeric ones, an optimal derivative step size can be derived (Heath 2002); the rule of thumb is a step size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$dx/x \sim \sqrt{\epsilon }$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon $$\end{document} is the relative accuracy of the function evalution. If this is similar to the machine precision, then in single precision \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$dx/x \sim 10^{-4}-10^{-3}$$\end{document} .

In many algorithms the rate equations are coupled to other non-linear equations (e.g. continuum radiative transfer equations), necessitating Newton–Raphson also for the case of optically thin lines.

Convergence criteria and issues

A universal convergence criterium is that the relative corrections should have become smaller than some fraction \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _n$$\end{document} ; \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\Delta n_{i}/n_i| < \alpha _n$$\end{document} . The SUMO code for example uses \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _n=0.01$$\end{document} (the largest relative population change between iterations must have fallen below 1%). In Monte Carlo simulations, statistical noise can sometimes make high-lying rarely interacting states never reach these kind of convergence levels. A simple handling of this is to limit the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _n$$\end{document} condition to a group of low-lying states exclusively. For example does Lucy (2003) apply the condition to the 5 lowest-lying states only, albeit with a more stringent value \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _n=10^{-5}$$\end{document} .

One may also pose a condition on the maximum allowed error, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$max{\left( \varvec{e}\right) }$$\end{document} , relative to the inflow and outflow rates (LHS and RHS of Eq. (19)). It sometimes happens that the errors have reached a small regime, but the NR stepping is jumping between two solutions bracing the true one. One may then consider a convergence criterium based on this second condition alone, possibly flagging the solution as having a larger error than desired. One can also attempt to change the NR step size to break the jumping loop.

The matrix equation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{R} \cdot \textbf{n} = \textbf{s}$$\end{document} (or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{J} \cdot \mathbf {\Delta n} = -\textbf{e}$$\end{document} ) has either a single solution, no solution, or an infinite number of solutions. The latter two cases occur if the matrix is singular, i.e. if at least one equation corresponds to another equation scaled with a (non-zero) constant. If that constant is non-unity, there is no solution. If it is unity, there is an infinity of solutions.

True matrix singularity, for real physical processes, is not expected to occur unless some term is forgotten or something is formulated incorrectly. However, limitations to the numeric precision of the computer can lead to two types of problems (e.g. Press et al. 1992): 1) “False" matrix singularity due to machine precision limitation. 2) Accumulated roundoff errors, giving a solution with unacceptably large errors. Both types of problems are more likely to occur the closer to singular the true matrix is. The second problem also gets more severe for larger equation systems, and large dynamic range (difference between largest and smallest values of solution).

In addition to this, even if the system is neither singular nor close to singular, issues with how blocks of the full equation system are split out and iterated with each other can lead to a (single) non-physical solution. The most obvious case of this is when the level population vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{n}$$\end{document} obtains one or more negative entries. Obtaining negative solutions is fundamentally a symptom of the same underlying problem as obtaining no solution—the equation system describes constraints that have no physical solution. Such a situation can arise if some quantities, external to the block, are unable to make the necessary adaptations to keep up with other changes in the system. It is therefore clearly important how one decouples the full problem into an iteration between blocks. Another way negative solutions can arise is when the linearization overshoots the solution, e.g. if the starting condition is too poor, or the derivative not accurate enough.

This situation is one motivation for co-solving the level populations of all the ions of a given element, instead of one by one. Such an approach is implemented by both Shingles et al. 2020 and van Baal et al. (2023)—here negative level populations are rare, whereas in SUMO, which solves ion by ion, they can happen quite often and corrective action become necessary. There is, however, always something that is held fixed whichever equation block one is solving. “No solution”, based on real issues with the equation system rather than numeric accuracy-related ones, can therefore always arise and needs to be dealt with in the codes. Examples of strategies are (1) stay with the previous block solution (and hope that the issues resolved in subsequent iterations), (2) reset to some specified solution (“reset initial guess”), and (3) reduce the step size for some quantities. One may also (4) attempt to backtrack on updates to other blocks. For example—instead of using the new photoionization and photoexcitation rates as estimated in the previous radiation field solution, one may try a step between the previous values, that gave a valid solution, and the new ones.

Because classes of algorithms come in a large number of flavours—and small changes can make a big difference, depending on application, there is basically an infinity of choices. Not only that—the performance of any given algorithm (stability, run time) varies over a practically infinite space (position-dependent density and composition, as well as starting guesses for all the parameters). Therefore, what is perhaps more meaningful than trying to compare and rank algorithms, is to devise a layered approach in the algorithm design, where a variety of fall-back strategies are resorted to in case the currently attempted one fails. At low level, this could for example be to retry with a different step size. At high level, it could be to redefine the whole equation block currently being attempted to solve for.

Coupling to the radiation field

The rate equations involve radiation field quantities; photoexcitation rates which depend on the line-averaged mean intensities \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{J}_{ij}$$\end{document} , i and j denoting the two levels involved, and photoionization rates which depend on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_\nu $$\end{document} convolved with the cross sections \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma _\nu $$\end{document} for the given bound-free transition. In regular Lambda iteration one completely decouples solutions to level populations and to radiative transfer, i.e. for each level population solution the radiation field terms are held constant. While this is a method with poor convergence properties for high-density situations and with resolved line transfer, the situation is much more benign in Sobolev modelling, especially for the low-density tail phase. This is because the Sobolev approximation treats multiple line scatterings as single ones, making the problem equation-wise simpler and numerically much easier to solve by iteration. It gives mostly robust and fast convergence for both SNe and KNe. The SUMO code implements a flexible way to select the degree of coupling by allowing the user to specify which (low-lying) levels are to be solved for in NLTE, and which (high-lying) are present for on-the-fly scattering and fluorescence but do not participate in collisionally induced emission. By reducing the NLTE coupling in this way can the Lambda iteration be used up to quite high densities, at least for Sobolev modelling.

For diffusion-phase modelling of SNe and KNe, better coupling than in classic Lambda iteration is preferable and sometimes necessary. If the radiation field is discretized into M frequencies, and the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_k$$\end{document} terms ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=1,M$$\end{document} ) are treated as variables to be solved for jointly with the level populations, the NR equation system (N rate equations plus M radiative transfer equations) takes the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{bmatrix} j_{1,n_1} & ...& j_{1,n_N}& j_{1,J_1}& .. & j_{1,J_M}\\ .. \\ j_{N,n_1} & ...& j_{N,n_N} & j_{N,J_1} & ..& j_{N,J_M} \\ t_{1,n_1} & ...& t_{1,n_N}& t_{1,J_1}& .. & t_{1,J_M}\\ .. \\ t_{M,n_1} & ...& t_{M,n_N} & t_{M,J_1} & ..& t_{M,J_M} \\ \end{bmatrix} \cdot \begin{bmatrix} \Delta n_1 \\ .. \\ \Delta n_N \\ \Delta J_1 \\ .. \\ \Delta J_M \\ \end{bmatrix} = \begin{bmatrix} -e_{n_1} \\ .. \\ -e_{n_N} \\ -e_{J_1} \\ .. \\ -e_{J_M} \\ \end{bmatrix} \end{aligned}$$\end{document}

where j denotes Jacobian entries from the rate equations and t correspondingly from the radiative transfer equations. One may also add \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_e$$\end{document} and T as variables by the two equations of charge and energy conservation (e.g. Hillier 1990). The issue is that the transfer equations depend on the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta n$$\end{document} population changes of all elements and levels, so for this coupling to give real convergence improvement it is preferable that the N rate equations include all levels for all species. This puts severe constraints on the number of levels per species one can treat.

In one of the pioneering works to address the matter-radiation coupling, Hillier (1990) used the transfer equation to eliminate the radiation field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta J$$\end{document} terms in the expressions involving the level population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta n$$\end{document} terms. In principle does the radiation field depend on population \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta n$$\end{document} terms at all depths. Retaining such a full dependency corresponds to “complete linearization” where one solves for population corrections at all depths simultaneously. It quickly becomes a prohibitive approach for any realistic SN/KN modelling problem, as one has tens of thousands of levels or more, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gtrsim $$\end{document} 100 zones even in 1D, giving an equation system exceeding size of a million and having \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$>10^{12}$$\end{document} entries. It can however be used for small, specific problems (e.g. Hillier et al. 1983; Hillier 1989).

A “local” operator limit corresponds to retaining dependencies of J in a given cell only on the level populations in that cell. Then the level populations can be solved cell-by-cell. This is attractive from the perspective of ease of parallelization, and also for a straightforward extension to multi-D. It is the operator used by CMFGEN for supernova applications. For non-local coupling, the simplest scheme is nearest-neighbor, which in gives a block-diagonal equation system. It is used by CMFGEN for stars (Hillier, priv. comm.).

Monte Carlo codes are by construction doing the radiative transfer in a separate computational step, and cannot direcly make use of the equation-based coupling schemes discussed above. However, they can achieve good convergence by other means. One key technique is the use of indestructible energy packets (Lucy 1999, 2005). The indestructability corresponds to enforcing radiative equilibrium which must hold at convergence, and enforcing this from the beginning improves the convergence situation. Another method to improve coupling in NLTE is the use of so called macro-atoms (Lucy 2002). An extensive review of these and other techniques in MC transfer is given by Noebauer and Sim (2019).

Photoexcitation rates.

Photoionization involves a continuum state (the ionized electron) and the rates therefore depend on integration of the radiation field over a quite broad range of frequencies. As such, it is relatively straightfoward to produce good estimators for photoionization rates with low or moderate resolution transfer equation solving, and in Monte Carlo simulations with small or moderate packet counts.

Photoexcitation, however, depends on the radiation flux over a very small frequency range; that of the line profile whose thermal width is of order 1 km/s (so \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta \nu /\nu _0 \sim 1/3\times 10^{-5}$$\end{document} ). Here, things get trickier. It is, for example, by no means obvious that a Monte Carlo simulation will have produced at least one packet interacting with each line in each cell. By not relying on counts but instead using information from all packets passing through lines (whether they interact or not) the statistics is improved (Lucy 1999)—but nevertheless one cannot be sure packets make such pass-throughs in all lines. Here, it is suitable to make checks whether rates have been generated at all, and if not perhaps choose to use the previous (non-zero) rate estimate instead. Or, if there are too many such instances, increase the packet count in the simulations. This illustrates a general aspect of Monte Carlo modelling that it comes with need for some special considerations compared with solving the transfer equation. On the other hand, solving the transfer equation without use of the Sobolev approximation requires resolving each line which puts much higher requirements on the number of frequency points.

A related issue is that storage of the photoexcitation rates for all transitions in each cell can become prohibititely expensive in 3D. Therefore, 3D modelling solving rate equations is so far limited to either ignoring these rates (Botyánszki and Kasen 2017; Botyánszki et al. 2018; van Baal et al. 2023) or estimating them from a low-resolution model of the radiation field (Shingles et al. 2020). Accurate treatment of photoexcitation rates in 3D NLTE modelling is a yet unsolved problem.

Further aspects

Starting guesses. In root finding with multi-dimensional Newton–Raphson it is important to have a good starting guess. Borrowing from Numerical Recipes: “Carefully crafted initial estimates reward you not only with reduced computational effort, but also with increased understanding and self-esteem”. This last point is valuable to regularly remind ourselves us, as the steadily increasing power of computers and algorithms may lead to loss of understanding of the basic physical behaviour of the system unless we stay involved.

For problems that are somewhat close to an LTE situation, LTE populations are the natural choice for the starting conditions. At low enough densities, however, LTE fails so severely that it is not useful for this purpose. Instead, expectations for the rough behaviour of the system, coming from a combination of theory and experience, provide the best starting guesses. For example, we know from previous solutions computed in the literature that nebular-phase CCSNe have ionization degrees of a few percent. We can combine this with the basic property that levels with forbidden deexcitation paths only will have much higher abundances than levels with at least one allowed deeexcitation path. Figure 7 illustrates this behaviour. A rough starting guess in the low density limit for a level i is then \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_i \approx \left( \sum _{j=1}^{i-1} A_{i,j}\right) ^{-1} \times max\left[ n_g C_{g,i}(T),\alpha _i n_e n_+\right] $$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _i$$\end{document} is the recombination rates, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_{g,i}$$\end{document} is the collisional excitation rate from ground (“g”) to level i, which scales with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\exp {\left( -kT/E_i\right) }$$\end{document} (see Jerkstrand 2017, for more details).

Specific software for linear equation system solving. The BLAS/LAPACK libraries provide robust and fast algorithms for basic linear algebra operations like solving a linear equation system. DGESV is the standard routine when the matrix has no particular simplifying properties to its structure (e.g., tridiagonal, symmetric). It employs LU-decomposition, which is a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^3$$\end{document} operation, and therefore becomes slow for large equation systems. However, as Fig. 8 shows, in practice the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^3$$\end{document} scaling will not become problematic in until N gets quite large, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \gtrsim $$\end{document} few hundred. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\lesssim $$\end{document} 500, the LU decomposition is subdominant to other data operations and the scaling is close to linear. For higher N, the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^3$$\end{document} scaling is closely reproduced. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\sim 500$$\end{document} is close to a sweet-spot for balance between accuracy and compute time. For most elements up to iron one is also close to the ionization edge at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \sim 500$$\end{document} , in full LS-coupling, so increasing this further is not expected to increase accuracy by much. At the same time, reducing N to lower values will not give significant total run-time improvements, especially when surrounding operations (loading the matrices) are considered as well. For trans-iron elements, especially open f-shell ones, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\sim 500$$\end{document} may not be sufficient and avoiding the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^3$$\end{document} regime requires application of techniques such as superlevels (more below).Fig. 8. Solution time for DGESV (red squares) and build times (blue stars) with randomly generated values in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{R}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{s}$$\end{document} , versus equation system size N. Also shown are some build times for different elements in a SUMO model (orange circles). Computed on a MacBook Pro Intel core i9 with gfortran compilation

There exists an advanced version of DGESV, called DGESVX, that may be a suitable second ansatz in a layered strategy. DGESVX allows for equation rescaling, iterative refinement, and other features that may enable solutions when DGESV fails for close-to-singular matrices.

Various operations and pre-conditions can improve the numeric situation. For example, CMFGEN performs two operations prior to the LU-decomposition solver step (Hiller, priv. comm.); 1) Recasting of the equations to solve for relative steps rather than absolute ones. 2) Pre-conditioning of the matrix with DGEEQ of LAPACK. Another example is that ARTIS solves for departure coefficients rather than absolute abundances. The value of such operations can be understood from inspection of the level populations example in Fig. 7: solutions can span 15 orders of magnitude in NLTE which can pose severe numeric challenges if the equations are solved for in “raw” form.

Iterative methods. An alternative to directly finding the exact solution to a linear equation system is to iterate approximate solutions until sufficient accuracy is reached. There is a vast array of methods and the most suitable one depends on the properties of the particular matrix. For rate equations, we have no general simplifying properties of the structure. In such cases, among the non-stationary methods, Conjugate Gradient on the Normal Equations (CGNE) or Generalized Minimal Residual (GMRES) are the most relevant to use.

So far, in no numeric approach to SN/KN spectral synthesis, that the author is aware of, have iterative methods been implemented. This is because N has been kept below \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 10^3$$\end{document} , and as Fig. 8 shows, in that regime 1) The solve time is still in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\propto N$$\end{document} regime not quite having transitioned to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^3$$\end{document} yet, and 2) Matrix build times are anyway comparable to the solve time. However, as computational capabilities improve, solutions for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\gtrsim 10^4$$\end{document} systems will become relevant, in particular for KN modelling, and for these LU decomposition may become replaced by iterative methods as the standard choice.

Superlevels. The number of levels, summed up over all elements and ion stages, quickly becomes large. This can give rise to challenges, the nature of which depends on the algorithm. For example, CMFGEN solves for all level populations in a cell simultaneously to allow for coupling between them (through the radiation field) (Hiller, priv. comm.). At full (fine-structure) level resolution, the matrix therefore becomes much too large to even store. SUMO, on the other hand, solves ion by ion. Storing is then not a problem, but nevertheless the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N^3$$\end{document} scaling of LU-decomposition makes large N models costly when N gets large. The situation becomes especially problematic for lanthanides which can have tens of thousands of levels all well below the ionization threshold.

One approach to address this is to use superlevels (Anderson 1989). This means that one bundles certain groups of levels into a single effective level. Typically one assumes an LTE distribution within the superlevel. Looking at the transitions between the fine-structure levels in a multiplet, typically the radiative rates are low and the collisional rates large. These conditions favor relative LTE distributions.

As an example of the scale reduction that can be achieved by superlevel use, the full-composition Type Ia models in Blondin et al. (2023) have 10,605 levels grouped into 2,338 superlevels—saving a factor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 16 in storage and up to a factor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 64 in matrix inversion time.

One should note that if one wants to retain individual, wavelength-separated lines connected to the superlevel, the energy discrepancies arising in the thermal balance (using the single specified energy of the superlevel) and radiative transfer (using several energies all slightly different from the single superlevel one) needs addressing. Hillier and Dessart (2012) discuss techniques for this for both hot wind and supernova applications.

Radioactive powering

The radioactive decays in SNe and KNe inject high-energy particles - gamma rays, electrons, positrons, and, for the heaviest r-process elements, also \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} -nuclei and fission fragments. These particles define the de-facto power source for the nebula already from an early point in time—the energy deposited by the explosion itself is quickly adiabatically degraded and is not able to create much radiation (the only exception to this are the Type IIP SNe in which the large radii of the red supergiant progenitors make explosion-deposited energy important for the first 2–3 months).

These high-energy particles transfer their kinetic energy to the ambient gas by the three processes of ionization, excitation, and heating. They may also lose energy by radiative losses such as Bremsstrahlung (this radiation may or may not be reabsorbed by the ambient gas), and by adiabatic losses. Ionizations by gamma-rays will eject bound electrons with quite high energies, called primaries. These can in turn cause further ionizations, these electrons are called secondaries. This cascading means a complete theoretical and numerical treatment of the physical situation is a significant challenge.

LTE codes do not solve rate equations and therefore information about what fractions of the energy loss that go into the various channels (ionization, excitation, and heating, and their specific subchannels) cannot be used. The only option in this case is to allocate all energy deposition to heating. If the degradation (in this context often called thermalization because all the energy is assumed to become heat) is fast and local, nothing needs to be solved for at all—the total radioactive decay power in each cell is simply used as a heating contribution term going into the temperature equation (Sect. 2).

Typically non-locality (transport) of gamma-rays are allowed for, however, as they can travel quite some distance in the ejecta, and in some cases partially escape already in the diffusion phase. Their main energy deposition mode is Compton scattering, in which typically a large fraction of the gamma-ray energy is transferred to an electron, creating a distribution of 0.1–1 MeV primary non-thermal electrons. The Compton scattering process is straightforward to solve the transfer equation for or to simulate with Monte Carlo methods.

Baryons and leptons, on the other hand, have much shorter mean-free-paths and are, in addition, quite easily trapped even by weak magnetic fields (Axelrod 1980). They are therefore often assumed locally absorbed.

As densities decrease with time as the nebula expands, at some point does the degradation start to be slow and/or non-local also for the leptons and baryons. Non-locality for leptons occurs only after about 5–10 years for Type II CCSNe, when NLTE is anyway necessary. An example is the situation in SN 1987A at an age of 8y, when positrons from ^44^Ti-decay power the metal-rich SN core and the spectrum probes their degree of spatial transport (Jerkstrand et al. 2011). For KNe, it was established in Barnes et al. (2016) that time-dependent effects (slow degradation) are important already after a few weeks or months.

LTE modelling

We will now study the different processes by which radioactive decay particles, and the particles generated by their interactions, deposit their energy to the ejecta. A sandbox to explore the various loss mechanisms have been coded up and is available at https://github.com/eliotayache/elec_degrad. It is useful to build intuition for how the different loss channels operate depending on energy and density, and which ones need consideration for a given problem at hand.

Ionization and excitation

The energy loss of a fast proton, or other nuclear projectile, due to collisional ionization and excitation of the ambient medium was worked out by Hans Bethe in the early 1930s, soon after the necessary quantum physics was in place. His work built upon the semi-classical foundation for the process established by Niels Bohr in 1913. The formula Bethe arrived at (see Fano 1963, for a detailed rederivation and discussion) is a good approximation for projectile energies much larger than the ionization potentials ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gtrsim $$\end{document} 500 eV, Axelrod 1980), and holds also for relativistic motion. This means it is useful for determining the distance (or equivalently, time) over which radioactive decay particles deposit the bulk of their energy (over 99.9% of the energy has been lost by the time a 1 MeV primary goes below 500 eV). Another property of the Bethe regime is that the energy loss occurs in a continuous manner, with small fractional losses for the primary in each ionization or excitation event. This property has consequences for analytic and numeric solution approaches.

With adaptation of the SI formula of Longair (2011) to cgs form,5 the Bethe formula is, for single species gas (the “np” superscript refers to nuclear projectile):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{dE}{dx}^{\rm{Bethe,np}}(E) = \frac{4 \pi q^4}{m_e} z_p^2 n_{eb}\frac{1}{v(E)^2} \left[ \ln {\left( \frac{2\gamma (E)^2 m_e v(E)^2}{\bar{I}}\right) }-\frac{v(E)^2}{c^2}\right] \ \text{ erg }\text{cm}^{-1}, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{eb}$$\end{document} is the number density of targets (all bound electrons per unit volume), q is the unit charge (which is 1 in cgs units, with dimension \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^{3/2}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {g}^{1/2}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document} ), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_p$$\end{document} is the (integer) charge of the projectile, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{I}$$\end{document} is a mean potential, considered as a parameter that typically needs experimental input for setting6. In 1933 Felix Bloch showed that a reasonable approximation is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{I} \approx 10\ \text{ eV } \times Z$$\end{document} , where Z is the atomic number. If that is inserted in the equation the formula is referred to as the Bethe–Bloch formula. The quantity dE/dx is sometimes called the stopping power. It is sometimes presented in units of erg  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {g}^{-1}$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^2$$\end{document} , obtained by division with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho $$\end{document} . The value of the constant factor is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4\pi q^4/m_e=7.34\times 10^{-10}\, \mathrm {erg\, cm^{4}\, s^{-2}}$$\end{document} . The corresponding energy loss per unit time is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L(E)\equiv dE/dt =dE/dx\times v$$\end{document} .

For electron projectiles, the formula is (again converting the SI form of Longair (2011) to cgs, and going to energy per time form):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} L^{\rm{Bethe,e}}(E) = \frac{2\pi q^4}{m_e} n_{eb} \frac{1}{v(E)}\left[ \ln {\left( \frac{\gamma (E) m_e v(E)^2 E_{\max }}{2\bar{I}^2}\right) } + f(\gamma (E))\right] \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(\gamma (E))$$\end{document} is a factor varying between 0.3 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-$$\end{document} 0.2 at the energies relevant here ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =1-3$$\end{document} ), and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{\max }$$\end{document} is the maximum transferable energy per collision, which for electron–electron collisions is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_{max} = \frac{\gamma ^2 m_e v^2}{1 + \gamma }$$\end{document} , so the loss can also be written

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} L^{\rm{Bethe,e}}(E) = \frac{4\pi q^4}{m_e} n_{eb} \frac{1}{v(E)}\left[ \ln {\left( \frac{E}{\bar{I}}\right) } -0.19 + g(\gamma (E)) \right] \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\gamma )$$\end{document} is a relativistic correction term only important for high energies.

Figure 9 shows the loss rate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^{\rm{Bethe,e}}(E)$$\end{document} for a target number density of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=10^7$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^{-3}$$\end{document} (the scaling is inverse with n). The loss rate has a minimum close to the characteristic injection energies of primaries, at or a bit below 1 MeV. The loss rate increases by a factor of a few going to lower energies, with maxima at about \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10\bar{I}$$\end{document} . We can see that for the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim $$\end{document} 1 MeV regime the relativistic correction term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\gamma )$$\end{document} may be ignored. Finally, note the moderate impact of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{I}$$\end{document} ; even a factor 10 change leads to a factor \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$<2$$\end{document} change in the loss rate, due to the logarithmic dependency.

For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} particles and fission fragments, we may use Eq. (24) instead of Eq. (26). These heavy particles are always in the non-relativistic regime.

Heating

Energy transfer to the free, thermal electrons occurs in a continuous manner as no bound states are involved - the long-range Coulomb potential leads to a smooth energy transfer process. Due to the repulsion and their small mass, the free thermal electrons are able to move away from the path of the projectile to some extent, which also contributes to this. As such we can assume their energy gain is not sufficient to make them leave the thermal pool and no cascading is obtained as in the ionization case.

The modelling of heating losses by fast particles has a long history in astrophysics, with the original motivation being the study of how cosmic rays impact nebulae. Spitzer and Scott (1969) derived a heating loss formula (for the secondary electrons created in the cosmic ray ionizations, see also Dalgarno et al. (1999):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} L_{e,\mathrm {heat}}(E)= \frac{4 \pi q^4}{m_e} n_{e,th} \frac{1}{v(E)} \Lambda (E)~\text{ erg}\;\text{s}^{-1}, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda (E)$$\end{document} is a dimensionless quantity tabulated in Spitzer (1962, a book, unavailable online). Note the similarity of the functional form to Eq. (24). A widely used formula for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda (E)$$\end{document} for SN/KN modelling is the result of Schunk and Hays (1971),

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Lambda (E)^{\rm{SH}} = \ln {\left( \frac{m_e v(E)^2}{\bar{h}\sqrt{4\pi q^2 n_{e,th}/m_e}}\right) }. \end{aligned}$$\end{document}

Note here the additional appearance of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{e,th}$$\end{document} . The Schunk and Hays (1971) formula corresponds to the non-relativistic limit of a more general expression derived by Inokuti et al. (1978),

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \Lambda (E)= \Lambda (E)^{\rm{SH}} + \frac{1}{2} \ln {\left( \frac{1}{1-\beta ^2}\right) } - \frac{1}{2}\beta ^2, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =v/c$$\end{document} . The Inokuti formula in turn builds from Fano (1963), but adds in collective excitations of the plasma. The additional terms in the (Inokuti et al. 1978) formula gives only small corrections below a few MeV (Fig. 9). Thus, for the purposes of SN/KN radioactivity modelling the Schunk and Hays (1971) formula should be sufficient.

The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda (E)$$\end{document} factor varies over the range \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$15-30$$\end{document} for relevant combinations of E ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 10^2-10^6$$\end{document}  eV) and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{e,th}$$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 10^3-10^{10}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^{-3}$$\end{document} ) for SN/KN applications. Thus, it is a reasonable approach to fix it to a constant, this is done e.g. in the KN modelling of Barnes et al. (2016).7

Discussion

Figure 9 plots the loss rates per target electron. For significantly ionized H and He-rich gases, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{e,th}$$\end{document} would be larger than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{eb}$$\end{document} , and heating is then the dominant loss mechanism for all energies. However, for heavy composition layers, models and observations indicate that the gas becomes at most a few times ionized (i.e., Fe would exist as Fe I-IV ions). In that case many more electrons remain bound than free (for example 23 bound and 3 free if all Fe is in the Fe IV state). The effective ratio of ionization/excitation to heating scales with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{eb}/n_{e,th}$$\end{document} (23/3 in the example) meaning most loss occurs by ionization/excitation in such a situation.Fig. 9. Electron energy loss rates as function of projectile energy. The collisional loss rates scale with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{\rm{target}}^{-1}$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{\rm{target}}=n_{e,th}$$\end{document} for heating and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{\rm{target}}=n_{eb}$$\end{document} for ionization/excitation, which is therefore embedded in the ordinate values for these. The adiabatic loss rate is just E/t (so does not depend on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_{targets}$$\end{document} )

Another point to take note of that the distribution of energy loss over the channels is not the same as the final distribution of energy going into the different channels. For example, if a 1 MeV electron transfers 100 eV to a bound electron in a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{I}=10$$\end{document}  eV atom, it has lost 100 eV by an ionization process, but the ionization energy of the gas has only increased by 10 eV. The other 90 eV resides now in another high-energy electron, and only a fraction of that will eventually be added on as further ionization energy as this particle cascades down.

Radiative losses

Radiative losses are for the most part small or neglegible compared to the collisional ones. The main radiative loss mechanism is Bremsstrahlung (Barnes et al. 2016). It can be described by (Seltzer and Berger 1986)

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} L_{\rm{rad,Brems}}(E) = r_e^2 \alpha \sum _i n_i Z_i^2 v(E)\left( E + m_e c^2\right) \phi (E,Z)~~\text{ erg}\;\text{s}^{-1}, \end{aligned}$$\end{document}

where summation is done over all ions i with charge \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z_i$$\end{document} . The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi (E,Z)$$\end{document} function is tabulated in Seltzer and Berger (1986) and has a typical value of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 10.

Formulation and numerics

Radioactive decay leads to injection of high-energy particles with both a total power and energy distribution being functions of time. For a given type of decay particle, this may be represented as a sum of decays per unit mass \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{n}_i(t)$$\end{document} at specific energies \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_i$$\end{document} , so total decay power \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\dot{\epsilon }(t)=\sum _i E_i \dot{n}_i(t)$$\end{document} erg \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {g}^{-1}$$\end{document} .

The energy loss rate per particle is given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{dE}{dt} = -L_{\rm{atom}}(E,t) - L_{\rm{elec}}(E,t) - L_{\rm{rad}}(E,t) -x(E)\frac{E}{t}, \end{aligned}$$\end{document}

where the last term is due to adiabatic losses (Waxman et al. 2018). Its pre-factor x has limits of 1 and 2 in the ultra-relativistic and non-relativistic limits, respectively. Adiabatic loss is not significant initially, but because it evolves as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^{-1}$$\end{document} compared to (for fixed gas state) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^{-3}$$\end{document} for the collisional loss terms, it will eventually come to dominate.

For a given decay particle of energy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_0$$\end{document} at creation, does a unique mapping \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{E}(E_0,t_0,t)$$\end{document} — the energy of the particle at any later time t have—exist? Such a mapping would come from solving the Eq. (31) IVP with initial condition \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E(t_0)=E_0$$\end{document} . However, the ionization, excitation, and heating losses depend not only on composition but also on the time-dependent gas ionization state, which in turn depends on the solution to the degradation evolution for all decay particles. Therefore, a full solution to this problem requires simultaneous time-evolving solution of the gas state (corresponds to determining the temperature in LTE) and the degradation equations. No such work yet exists in the literature. And clearly, the answer to the posed question is no—the answer is density and composition dependent, and in addition is the evolution dependent on all decay particles.

Anything of generic nature must therefore assume a fixed gas state, or a parameterized time evolution of the gas state. In the approaches in the KN literature so far those assumptions are

1. Simplified loss functions taken as constants independent of energy and time (Waxman et al. 2018).
2. All elements are in the singly ionized state (Barnes et al. 2016; Hotokezaka and Nakar 2020b).
3. All elements are in the neutral state (Kasen and Barnes 2019). The \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{E}_i(E_0,t_0,t)$$\end{document} functions can then be uniquely determined and the power at time t becomes, ignoring secondary electrons,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} P(t) = \sum _i \int _0^t \dot{n}_i(t')\left[ L_{i,\mathrm {atom}}(t',t) + L_{i, \mathrm {elec}}(t',t)\right] dt'. \end{aligned}$$\end{document}

The solution of Eq. (31), needed to compute Eq. (32), involves solving a first-order IVP with non-constant coefficients. Analytic limits based on the simplifications outlined above are extensively discussed in Barnes et al. (2016), Waxman et al. (2018, 2019), Kasen and Barnes (2019), Hotokezaka and Nakar (2020b). A public code to solve the equations for KN ejecta (Hotokezaka and Nakar 2020b) is available (Hotokezaka and Nakar 2020a).

Examples

Figure 10 (left) shows a calculation from Barnes et al. (2016) of the so called thermalization efficiency in a KN, defined as the power at time t (Eq. (32)) divided by the radioactive decay power at the same time. Observations of AT2017gfo covered up to 20 d in optical/NIR and yet later (up to 74 d) with the Spitzer mid-infrared telescope. Over observable epochs, thermalization is clearly a time-dependent phenomenon requiring solutions of Eqs. (31) and (32).

The right panel of Fig. 10 shows different models, with time-dependent thermalization, compared to observations of AT2017gfo, from Hotokezaka and Nakar (2020b). It shows that modelling the powering allows constraints to be put on the KN composition. Fig. 10Left: Thermalization efficiency versus time for a homologous expansion trajectory corresponding to typical KN densities ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho (1d)=8\times 10^{-15}$$\end{document} g \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^{-3}$$\end{document} ). Right: Observations of AT2017gfo compared to models with different nuclear compositions, and with time-dependent thermalization. The figure demonstrates that the late-time KN data is able to constrain the composition. Images reproduced with permission from [left] (Barnes et al. 2016) and [right] (Hotokezaka and Nakar 2020b), copyright by AAS

NLTE modelling: Boltzmann solutions

For NLTE modelling the specific power going into the different channels of ionization, heating and excitation must be solved for. The ionization channel breaks down into subchannels for each ion-electron shell pairs (e.g. O I K, O I L, O I M, O II K...), whereas the excitation channel breaks down into the set of considered bound-bound transitions for each atom/ion. Electrons at each energy in the degradation cascade will contribute differently to these different channels—the degradation energy spectrum therefore needs to be solved for. Free electrons in the ejecta span energies from the thermal population ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim kT \sim $$\end{document} eV) to the primary electrons released in radioactive decays or other high-energy process ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document}  MeV). Given the six orders of magnitude in the energy cascade, this is a challenging task.

Current codes solve the Boltzmann equation under the two assumptions

1. Steady state (no time-dependent degradation).
2. Locality (no spatial transport). This ansatz is valid if the particles lose their energy over temporal and spatial scales much smaller than those of changing physical conditions. The first such solutions in the SN context were by Meyerott (1978) and Axelrod (1980). Later works include those of Kozma and Fransson (1992); Xu and McCray (1991); Lucy (1991). The first full solutions for KNe were done with the SUMO code.

Whereas LTE modelling can bundle non-thermal energy transfer to excitation and ionization into a single treatment, in NLTE it needs separation. The use of the Bethe formalism with an empirical \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{I}$$\end{document} is therefore less useful, because those empirically determined values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{I}$$\end{document} involve the combined effects of excitation and ionization. In addition, even though the low-energy regime where the Bethe formalism starts to fail ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim $$\end{document} 500 eV) is unimportant for thermalization time/distance, it is important for the ionization and excitation rates needed in NLTE. For these reasons, one now needs cross sections for individual processes, also extending down to low energies ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 1–500 eV).

Ionization terms

Let us begin with the total cross sections for ionization (integrated over all energy transfer values). A commonly used semi-empirical formalism is that of Lotz (1967). In the initial work of Axelrod (1980), the following flavour of this formalism was used:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sigma ^{\rm{tot}}(E) = \frac{{\rm{const}}}{\beta ^2 m_e c^2}\sum _{i=1}^{N {\rm{shells}}} \frac{n_i}{I_i}\left[ \ln {\left( \frac{\beta ^2m_e c^2}{2I_i}\right) } - \ln {\left( 1-\beta ^2\right) }-\beta ^2\right] , \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {\text const}=6.9 \times 10^{-38}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {erg}^2$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^2$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_i$$\end{document} is the number of electrons in shell i per atom, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =v/c$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_i$$\end{document} is the ionization potential of shell i. Later works for SN modelling have used dedicated theoretical or experimental determinations for the total cross sections when available, but the Lotz (1967) formalism still finds use for KN modelling (Pognan et al. 2022b).

A complication arises for inner shell ionizations, because refilling of the created "hole" leads to Auger electron ejection and/or X-ray fluorescence. Different codes take different approaches to dealing with this. The SUMO code, for example, lets the inner shells contribute to the ionization cross section, and the excess energy is used to raise the ionization rate (which corresponds to treating all electrons like valence ones). In Shingles et al. (2020), an improved treatment is devised, implemented into the ARTIS code, by coupling in a feedback loop for the Auger electrons. No large impact on Type Ia model spectra were seen, likely because for iron-group elements the contribution by the inner K and L shells is not so large, with 10 electrons in total compared to the 16 M-shell electrons. The impact may be larger for both lighter and heavier elements, when the valence count is not so dominant as in the iron-group.

Differential cross sections. The ionizations lead to the creation of secondary high-energy electrons. To solve for the Boltzmann distribution, we need to know the energy distribution of these, i.e. not just the total ionization cross sections \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma ^{\rm{tot}}(E)$$\end{document} discussed above, but also the differential ones \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma (E,\epsilon )$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon $$\end{document} is the energy transfer. Note the unit for these is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^2$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {erg}^{-1}$$\end{document} .

Stopping power depends on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\int \sigma (E,\epsilon ) \epsilon d\epsilon $$\end{document} , i.e. the first moment of the differential cross section. The Bethe formula for stopping power, Eq. (26) (valid for energy transfers \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I \ll \epsilon \ll E_0$$\end{document} ) derives from the first moment of a differential cross section following a "Rutherford" form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma \propto 1/\epsilon ^2$$\end{document} . This gives the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ln (E)$$\end{document} factor for the stopping power, because we need to integrate the loss over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1/\epsilon $$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1/\epsilon ^2 \times \epsilon $$\end{document} ).

More detailed information can be obtained from experiment. Opal et al. (1971) measured the secondary electron energies (up to 200 eV) resulting from the impact of high-energy electrons (up to 2 keV) on several molecules and inert atoms. The differential cross section was found to be well described by the formula

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sigma (E,\epsilon ) = {\mathrm {const}}(E) \times \frac{1}{1+\left( \epsilon /\bar{E}\right) ^2}, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{E}$$\end{document} was found from the experiments to lie in the range \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( 0.5-1\right) I$$\end{document} depending on element and projectile energy. This equation takes the Rutherford form in the high energy limit, but levels out in the low-energy limit. One should note that also inner shells may contribute electrons to this experimental ejection process, and also that in the experiment no distinction can be made between the two post-collisions free electrons, necessitating a definition of secondary electrons as the less energetic post-collision one.

The primaries in the Opal experiment had much lower energies ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lesssim $$\end{document} 200 eV) than the primaries in radioactive decays ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 1 MeV). For high-energy primaries, Manson et al. (1975) reported differential cross sections for proton impact on He, up to the MeV range. The differential cross section is found to have a distinct steepening for high secondary energies compared to Eq. (34); for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_0=100$$\end{document} keV impact that steepening begins at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon \approx 100$$\end{document}  eV, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_0=300$$\end{document} keV at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon \approx 300$$\end{document}  eV, and for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E_0=1$$\end{document}  MeV at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon \approx 1$$\end{document} keV; thus the steepening occurs uniformly at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon /E_0 \approx 10^{-3}$$\end{document} . A power law fit to the steeper part has an exponent of between − 4 to − 5. Assuming the behaviour is similar for electron projectiles, one has a prescription available for the differential cross section over the whole energy interval for the projectile.

Excitation terms

In excitations the energy transfer is fixed to the transition energy \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta E$$\end{document} . What is needed is just the cross section for the process as function of projectile energy. For large projectile energies this is, in Bethe theory for electron collisions, given by (Mott and Massey 1949; van Regemorter 1962)

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sigma ^{tot}(E) = 3 \pi q^4 f_{osc} \frac{1}{E \Delta E} \ln {\left( \frac{4E}{\Delta E}\right)}, E \gg \Delta E \end{aligned}$$\end{document}

This expression is used for all energies in the Boltzmann solver developed by Kozma and Fransson (1992) and used by the SUMO code, for those transitions lacking specific cross section data. Note that the Bethe limit is more readily obtained for line transitions than ionization continua as the line transitions have lower energy.

For light elements (up to iron) the energy going into excitations is typically comparable to, or smaller by a factor few, than the energy going to ionizations (Kozma and Fransson 1992). At the relatively neutral conditions in nebular CCSNe, which leads to the thermal electrons taking less than \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim $$\end{document} 50% of the power, the emissivity from non-thermal excitation can therefore be significant.

For r-process elements, no calculations yet exist. However, one may comment that the ionization state of KN ejecta appears significantly higher than in CCSNe (Hotokezaka et al. 2021; Pognan et al. 2022b). This means that direct emission from non-thermal ionizations and excitations carries less power compared to cooling emission following the non-thermal heating (Sect. 2). While the non-thermal ionizations are crucial for setting the ionization balance and determining which ions emit the cooling, non-thermal excitations have no such indirect impact. On the other hand, some of the r-process elements have orders of magnitude more line transitions than any light elements — which would boost the importance of non-thermal excitations. A clear picture of the overall effect awaits detailed calculations.

Solution approach

(Lack of) thermal electron coupling. If the thermal electrons ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$E \lesssim $$\end{document} 1 eV) are to be included in the Boltzmann solution, one must allow transfers that both decrease and increase the primary electron’s energy — and also account for the creation and removal of electrons by ionization and recombination terms. There will be coupling between all energies in general, i.e. a “full” matrix. There are no published solutions for such an approach in the literature.

Instead, all published solutions treat the thermal electrons separately. This is possible because quite few of the secondary electrons have such low energies that they are in the thermal regime at creation. For the Opal prescription of secondary electron distributions (Eq. (34)) and a standard choice \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{E} = 0.6 I$$\end{document} , the low-end 10th percentile of the distribution is at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 0.2I$$\end{document} . Then even the lowest ionization potentials ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I \sim 6$$\end{document}  eV) give less than 10% of secondaries below 1 eV. We should therefore be able to put a clear distinction, around 1 eV for typical conditions, between a thermal electron group below this threshold and a non-thermal one above.

Formulation and numerics. Let z be the flux of non-thermal electrons at energy E per unit area and energy ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^{-2}$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {erg}^{-1}$$\end{document} ). It is related to the spectral number density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_E$$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^{-3}$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {erg}^{-1}$$\end{document} ) by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z = n_E v(E)$$\end{document} . If the source term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$S_E$$\end{document} (containing injection by radioactive decay) is given in units of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^{-3}$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {erg}^{-1}$$\end{document}   \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {s}^{-1}$$\end{document} , the steady-state distribution for z is determined by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&S_E(E) + \sum _i n_i \left[ \int _{I_i}^{E_{\max }} z(E') \frac{ d \sigma _i^{\mathrm{ion+exc}}}{d \epsilon }(E',E'-E) dE' \right. \nonumber \\&\left. + \int _{E+I_i}^{E_{\max }} z(E') \frac{ d \sigma _i^{\rm{ion}}}{d \epsilon }(E',I_i+E) dE'\right] \nonumber \\&= z(E) \sum _i n_i \int _0^E \frac{ d\sigma _i^{\mathrm{ion+exc}}}{d \epsilon }(E,\epsilon ) d \epsilon - \frac{d}{dE}\left[ z(E)\frac{dE}{dx}_{\rm{heat}}(E)\right] \end{aligned}$$\end{document}

Here, the first term on the LHS represents injection by radioactive decays, the second term injection into the bin due to downscattered primaries, and the third term creation of secondaries following ionization events. On the RHS, the first term is loss from the bin due to ionizations and excitations, followed by a term representing heating losses (see Xu and McCray 1991, for derivation). This last term, not treated by Spencer and Fano (1954), changes Eq. (36) from an integral equation into an integro-differential one. Note that for excitations, the differential cross section can be treated as a delta function in the integral, in practice creating a separate sum term for these events.

Radiation terms. Including radiative losses, such as Bremsstrahlung or synchrotron, can be done by adding further (negative) source terms, which will be of the same general form as the heating loss term as long as the radiation process can be described as continuous. Self-absorption of this radiation may conceivably be added in if it creates non-thermal electrons. Neither radiation process is expected to be important for SNe, whereas for KNe Bremsstrahlung may play some role (Barnes et al. 2016). To date, radiation terms have not been implemented in any solutions in the literature.

Solution approach. As long as the thermal electron pool is uncoupled, and energy gain by the non-thermal electrons (by collisions with higher-energy ones) is ignored, existing electrons can only move downward in energy. Eq. (36) can therefore be solved from injection energy down to the lowest treated energies in single steps, not requiring matrix storage and inversion. It is a first-order integro-differential IVP, solvable by e.g. implicit Euler.

The discontinuous nature of bound-bound and bound-free transitions (e.g. the excitation cross section in Eq. (35) discontinuously drops to zero at the threshold energy as the peak lies below threshold at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.68 \Delta E$$\end{document} ) calls for some special care in the numerics. One should also pay attention to the fact that the differential ionization cross sections are only large at transfer energies of order of the ionization potential. Therefore, the integrations in Eq. (36) can involve integration from e.g. 10 eV to 1 MeV, using a kernel function that is only large up to say 100 eV for the energy transfer. Resolving this integral over the energy transfer argument means that the primary energy must be resolved to the same precision. In short: it will be problematic to use anything but a linear energy grid, and the resolution of this must be better than the ionization potentials.

It turns out that if the output needed are the various fractions of energy going into different channels, that can be calculated by an alternative formalism where one steps up in energy instead of downwards (Xu and McCray 1991; Lucy 1991). In this formalism the degradation spectrum z(E) is never explicitly solved for. One advantage of this approach is that the fractions will converge towards final values quasi-independent of the injection energy. Thus, there is a possibility to stop the integration with some criterion on the change rates. It is also more efficient if one needs to compute the solution many times, for different injection energies. However, different types of numeric issues show up in this approach compared to the degradation solution approach—see Sect. 4 in Xu and McCray (1991) for an in-depth discussion.

Some results

The left panel of Fig. 11 shows the level population boosts in the metastable states of He I (emitting the important 1.083 and 2.056 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} m lines) when non-thermal excitations and ionizations are considered (Lucy 1991). The departure coefficients reach \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log ~b = 4-5$$\end{document} which illustrates how powerful these mechanisms can be for certain elements and levels.

The right panel shows the fraction of radioactive power going into different channels as function of ionization state of the gas. We see how heating goes over 50% already at electron fractions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_e \equiv n_e/n_{\rm{atoms}}\gtrsim $$\end{document} 3%. This is why temperature calculations are not very sensitive to the details of non-thermal degradation in the diffusion phase ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_e \gg $$\end{document} 3%) or even in the nebular phase of CCSNe ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_e \sim $$\end{document} 10%). We also see that the atomic stopping contribution is about equally distributed between ionization and excitation for iron.Fig. 11Left: Departure coefficients for He in a Type Ib CCSN model at 15d (photospheric phase) when non-thermal effects are considered. Right: Fraction of non-thermal energy going into channels of heating, ionization and excitation, versus electron fraction in a pure Fe I + Fe II gas (so \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_{\rm{FeII}}=x_e$$\end{document} ). Images reproduced with permission from [left] (Lucy 1991) and [right] (Kozma and Fransson 1992), copyright by AAS

Summary and outlook

We have reviewed some of the key aspects for computing spectral models of explosive transients; supernovae and kilonovae. The emphasis has been on the computational techniques and considerations, and to sort out which approximations are in use and what they mean.

The main goal of spectral modelling is to, eventually, be able to infer the composition of SN and KN ejecta from observations. Figure 12 shows an overview of where the field stands in its current ability to do this.Fig. 12. Diagnostic situation for elements in the periodic table from SN and KN spectral modelling as applied to observational dataImage reproduced with permission

For the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Z=1-30$$\end{document} elements, SN spectral studies since the early 1980s have by now brought about a relatively mature picture of which signatures we observe and do not observe in spectra, and what we can infer from them. The elements are divided into the three categories of “good diagnostic potential”, “moderate diagnostic potential” and “challenging to diagnose”. The first category typically has at least two clear lines whose formation is relatively well understood and is abundance constraining. The last has either no clear lines, or perhaps a single line that is however sensitive to model uncertainties. The middle category is somewhere in between. It may be seen as a success of the field over its first 45 years that 18 out of these 30 elements are now blue or green (for recent efforts see e.g. Dessart et al. 2021; Fang et al. 2022; Liljegren et al. 2023; Barmentloo et al. 2024; van Baal et al. 2024). This accomplishment has been made possible by the combination of high-quality spectral data and the development and application of synthesis codes addressing the various challenges outlined in this review.

For the trans-iron group elements, spectral searches and analysis began only with the first KN in 2017. It is at this point not possible to meaningfully assess the diagnostic potential of many of these elements—we have neither sufficient observations nor accurate enough spectral models with r-process elements either in SNe or KNe. However, for KNe this is now rapidly coming about. I have marked in red the r-process elements that in the current literature have been identified as plausible candidates for spectral features in AT2017gfo, or appearing to have potential for identification in KN spectra (Watson et al. 2019; Domoto et al. 2021, 2022; Hotokezaka et al. 2022, 2023; Sneppen and Watson 2023; Pognan et al. 2023, 2025; Gillanders et al. 2024; Banerjee et al. 2025; Vieira et al. 2023).

To increase the blue and green in this diagram, on the radiative transfer modelling side we need to steadily improve on the physical treatments of the various foundational processes — hopefully this review contributes to visualizing the roadmap for this endeavor. Equally important is the development of accurate atomic data that goes as input to these model (in particular for the trans-iron group elements), development of 3D hydrodynamic models for various types of explosions and mechanisms, and high-quality observations over the whole UV to mid-infrared range and for the different epochs probing different layers of the ejecta. By the combined and linked efforts in these categories we can be optimistic to eventually be able to determine the origin of elements as directly observed at their production sites in supernovae and kilonovae.