© S. Blondin et al. 2022

1 Introduction

Accurate radiative-transfer (RT) calculations remain a key challenge in the study of astronomical transients such as supernovae (SNe). While advances in computational capabilities and theoretical understanding have allowed great progress in the ability to simulate radiation transport, the large number of physical processes involved, in particular opacity from thousands of atomic transitions with a mixed absorptive and scattering character, prohibit comprehensive 3D calculations based on first principles. Several physical approximations of different forms , in particular different treatments of the significant deviations from local thermodynamic equilibrium (LTE) , are employed by different RT codes to calculate the properties of the gas and of the radiation field. Approximate treatment of atomic physics is also required due to the partially calibrated atomic data.

The back-reaction of radiation on the hydrodynamics provides an additional challenge, requiring the solution of hydrodynamic equations coupled to the RT solution. However, in many cases, and in particular for Type Ia SNe (SNe Ia) at phases beyond several days, which we focus on here, the radiation carries a negligible fraction of the energy and the ejecta are freely expanding homologously. The RT problem in these cases is decoupled from the hydrodynamics problem, the latter providing the initial ejecta profiles (‘ejecta models’ hereafter). The ejecta profiles include density, composition, and initial temperature as a function of position. The initial time (of order 1 day) is such that on the one hand it is much larger than the explosion timescale (of order 1 second) so that the expansion is nearly homologous and on the other hand sufficiently early such that radiation has hardly diffused across the ejecta and the only evolution is due to adiabatic expansion and radioactive decay.

Comparisons of the results of different RT codes for the same ejecta models play an important role in estimating the accuracy of different approximations and can be used to validate new codes. The number and sophistication of RT codes has significantly developed in recent years, increasing the need for diverse benchmark ejecta models that will allow detailed and careful comparisons. In this paper, we describe the first results of a collaborative effort of ten groups around the world that are developing existing RT codes (in alphabetical order: ARTIS, CMFGEN’ CRAB, KEPLER, SEDONA, STELLA, SUMO, SuperNu, TARDIS, and URILIGHT; see Sect. 3 for descriptions and references) to create a systematic code-comparison framework. As a first important step, all groups agreed on a set of four test model ejecta and standardised output formats. Each group calculated the resulting radiative display with their respective codes for the same ejecta models and shared the results in a new public electronic repository on GitHub¹. We did not attempt to agree on a specific setup for each code, nor did we synchronise atomic data between codes.

The structure of the paper is as follows. In Sect. 2, the benchmark models are described. These include two simplistic sub-Chandrasekhar-mass (sub-M_Ch) toy models with profiles that are defined analytically and two more realistic models that result from hydrodynamical simulations of the M_Ch delayed-detonation scenario. All the models considered here are spherically symmetric (1D). We give short descriptions of the RT codes that were employed in this first comparison in Sect. 3, with emphasis on the main physical approximations that are used in each. In Sect. 4, we describe the publicly available repository of results. In particular, descriptions are provided of the output files and of Python codes that are included for reading them. In Sect. 5, several example comparisons of the results of the different codes are provided in order to illustrate the available outputs. We voluntarily make no attempt to analyse the sources of discrepancies. While the comparisons focus on observable aspects of the emission, comparisons to observations and conclusions regarding the implications for the applicability of the codes are intentionally not addressed in order to keep the focus of the paper on the description of the comparison. Finally, we outline future plans for this comparison project in Sect. 6.

2 Test models

The code-comparison test suite consists of four SN Ia models. Two are sub-M_Ch models with analytic density and composition profiles (Sect. 2.1; ‘toy’ models), and the remaining two are M_Ch models resulting from hydrodynamical simulations (Sect. 2.2; DDC models). The models were set up or selected based on their ⁵⁶Ni yield, in order to have two models corresponding to ‘normal’ SNe Ia (toy06 and DDC10 with ~0.6 M_⊙ of ⁵⁶Ni) and two low-luminosity models (toy01 and DDC25 with ~0.1 M_⊙ of ⁵⁶Ni). We present the toy and DDC models in turn in the following sections, and summarise their properties in Table 1. The density profiles at a reference time of 1 day post explosion and initial composition profiles are shown in Figs. 1 and 2, respectively.

2.1 Toy models

The toy models were generated using the mk_snia_toy_model.py script (see Sect. 4) using the --highni (for the normal SN Ia model) and --lowni (for the low-luminosity SN Ia model) options. Both models have a total mass M_ej = 1.0 M_⊙, a kinetic energy E_kin = 10⁵¹ erg, and are calculated at the time t_f = 2 days post explosion. The models have an exponential density profile (e.g. Jeffery 1999),

$$\rho ={\rho }_{\text{c}}{\text{e}}^{-v/v_{\text{e}}},$$(1)

where

$$v_{\text{e}}=\sqrt{\frac{E_{\text{kin}}}{6M_{\text{ej}}}}\approx 2895\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}$$(2)

is the e-folding velocity, and

$${\rho }_{\text{c}}=\frac{M_{\text{ej}}}{8\pi v_{\text{e}}^3{t}_{\text{f}}^3}\approx 6.32\times {10}^{-10}\text{g\hspace{0.17em}}{\text{cm}}^{-3}$$(3)

is the central density at the chosen time.

The ejecta consist of n spherical shells indexed by i = 1,…, n. We use a uniform velocity grid with width Δv = 50 km s⁻¹. Each shell then has an inner and outer velocity coordinate given by: v_0,i = (i − 1)Δv and v_1,i = v_0,i + Δv.

As in Jeffery (1999), we define the dimensionless radial coordinate z = v/v_e which we use to compute the mass of each shell as:

$$M_i={\displaystyle \int \rho \left(z\right)\text{d}V}={\rho }_{\text{c}}{\displaystyle \int e^{-z}\text{d}V,}$$(4)

where the integral runs from z_0,i to z_1,i and the volume element is $\text{d}V=4\pi r^2\text{d}r=4\pi v_{\text{e}}^3{t}_{\text{f}}^3{z}^2\text{d}z$, where we assume a homologously expanding ejecta (r = vt = v_ezt).

Likewise, each shell volume V_i is computed from the inner and outer radii (r_i,{0,1} = v_i,{0,1}t_f), which results in the mean density of each shell:

$${\rho }_i=\frac{M_i}{V_i}=3{\rho }_c\frac{e^{-z_{0,i}}\left(z_{0,i}^2+2z_{0,i}+2\right)-e^{-z_{1,i}}\left(z_{1,i}^2+2z_{1,i}+2\right)}{z_{1,i}^3-z_{0,i}^3}.$$(5)

Starting from the central shell, we keep adding successive shells until the total mass is within 0.01% of the required ejecta mass. For the 1.0 M_⊙ toy models considered here, this results in 807 shells, where the last shell has an outer velocity of 40 350 km s⁻¹. We note that not all codes use this velocity grid; in particular, the number of shells and maximum velocity vary significantly among codes (see Sect. 3).

We assume the ejecta consist of a central region composed of iron-group elements (IGEs), here only ⁵⁶Ni and its decay products ⁵⁶Co and ⁵⁶Fe, and an outer region composed of the intermediate-mass elements (IMEs) Ca, S, and Si with constant mass fractions throughout the layer: X(Ca) = 0.1, X(S) = 0.35, and X(Si) = 0.55. These values were chosen based on the delayed-detonation model DDC10 (see Sect. 2.2 below). Our ‘normal’ toy model consists of 0.6 M_⊙ of ⁵⁶Ni and 0.4 M_⊙ of IMEs, while our low-luminosity toy model consists of 0.1 M_⊙ of ⁵⁶Ni and 0.9 M_⊙ of IMEs. The ⁵⁶Ni and IME composition profiles are smoothly connected using an analytic function (here a cosine bell) over a mass interval ∆M_trans (set to 0.2 M_⊙ for the low-luminosity model and 0.4 M_⊙ for the regular model). At a given mass coordinate m, the ⁵⁶Ni mass fraction is set to:

$$X_{{}^{56}\text{Ni}}\left(m\right)=\{\begin{array}{ll}1\hfill & \text{for\hspace{0.17em}}m<m_1,\hfill \\ 1-f_{\cos }\left(m\right)\hfill & \text{for\hspace{0.17em}}{m}_1\le m\le m_2,\hfill \\ 0\hfill & \text{for\hspace{0.17em}}m>m_2,\hfill \end{array}$$(6)

where m₁ = M(⁵⁶Ni) − ∆M_trans/2, m₂ = M(⁵⁶Ni) + ∆M_trans/2, and

$$f_{\cos }\left(m\right)=\frac{1}{2}\left\{1-\cos \left[\left(\frac{m-m_2}{\text{Δ}{M}_{\text{trans}}}\right)\pi \right]\right\}.$$(7)

The IME mass fraction is then simply set to $X_{\text{IME}}\left(m\right)=1-X_{{}^{56}\text{Ni}}\left(m\right)$. Our toy models therefore consist of only six chemical species or isotopes (⁵⁶Ni, ⁵⁶Co, ⁵⁶Fe, ⁴⁰Ca, ³²S, ²⁸Si): this was a voluntary choice in order to keep the model as simple as possible while retaining the defining characteristics of a SN Ia.

The initial temperature profile at t_f = 2 days is determined by solving the first law of thermodynamics assuming a radiationdominated gas, local energy deposition from ⁵⁶Ni decay, and no diffusion (i.e. the temperature in each zone is solved independently of the adjacent zones). Given these assumptions, the temperature at t_f can be determined analytically by noting that the time-weighted internal energy, tE(t), is equal to the time-integrated time-weighted decay energy deposition rate, ∫tQ(t)dt, as noted by Katz et al. (2013):

$$T_i={\left(\frac{{\rho }_i{X}_i{\left({}^{56}\text{Ni}\right)}_0{\displaystyle {\int }_0^{t_{\text{f}}}t{q}_{\text{Ni}}\left(t\right)\text{d}t}}{a{t}_{\text{f}}}\right)}^{1/4},$$(8)

where X_i(⁵⁶Ni)₀ is the ⁵⁶Ni mass fraction at t ≈ 0 in the i’th cell, a is the radiation constant, and q_Ni(t) is the energy release rate per unit mass (ignoring neutrinos) of the ⁵⁶Ni→⁵⁶Co→⁵⁶Fe decay chain. In this formulation, we ignore the time-weighted internal energy shortly after explosion, E(t₀)t₀.

It is clear from Eq. (8) that the temperature is predicted to be zero in zones devoid of ⁵⁶Ni (≳12 000km s⁻¹). We therefore impose a constant temperature floor in these zones of 5000 K. The resulting excess internal energy is quickly radiated away because these layers have a relatively low optical depth, such that the impact on the longer-term radiative display is negligible.

Fig. 1 Density profiles of our model set at an adopted time of 1 day post explosion.

Fig. 2 Composition profiles at the start of our simulations in our model set (2 days post explosion for the toy models and ~1 day post explosion for the DDC models). For the toy models, this represents the full set of species present (56Ni and decay products 56Co and 56Fe, as well as IMEs: Ca S, Si), while for the DDC models only a subset of species are shown for illustration. The 56Ni mass fraction is given at the time of explosion. We also show the total IGE mass fraction (from Sc to Ni; dashed line) and the total IME mass fraction (from Na to Ca; dotted line). The total IGE mass fraction coincides with the 56Ni0 line in the toy models and is not shown for sake of clarity.

2.2 Delayed-detonation models

In addition to the two toy models above, we consider two M_Ch delayed-detonation models resulting from hydrodynamical simulations. These were chosen from the DDC model grid presented in Blondin et al. (2013) in order to closely match the ⁵⁶Ni yields of the toy models: Our low-luminosity model DDC25 yields ~0.12 M_⊙ of ⁵⁶Ni (cf. 0.1 M_⊙ for the toy01 model), and our ‘normal’ DDC10 model yields ~0.52 M_⊙ of ⁵⁶Ni (cf. 0.6 M_⊙ for the toy06 model). We refer the reader to Blondin et al. (2013) for a detailed description of the DDC models.

The outputs of the hydrodynamical modelling correspond to 30–60 s post explosion, by which time the ejecta have reached a state of homologous expansion. We applied a small amount of radial mixing to the hydrodynamical input with a characteristic velocity width ∆v_mix = 400 km s⁻¹ to smooth sharp variations in composition. The ejecta were then evolved to 0.5 day post explosion by solving the energy equation given by the first law of thermodynamics, assuming the plasma is radiation dominated and neglecting diffusion. Apart from the two-step ⁵⁶Ni→⁵⁶Co→⁵⁶Fe decay chain, we also treat eight additional two-step decay chains associated with ³⁷K, ⁴⁴Ti, ⁴⁸Cr, ⁴⁹Cr,⁵¹ Mn, ⁵²Fe, ⁵⁵Co, and ⁵⁷Ni, and a further six one-step decay chains associated with ⁴¹Ar, ⁴²K, ⁴³K, ⁴³Sc, ⁴⁷Sc, and ⁶¹Co (see Dessart et al. 2014).

The ejecta at 0.5 day are then remapped onto the 1D, non-LTE, radiative-transfer code CMFGEN of Hillier & Dessart (2012) and evolved until ~1 day post explosion (0.976 day for the DDC10 model and 1.3 days for the DDC25 model; see Sect. 3.2 for details), at which point the ejecta serve as initial conditions for the other radiative-transfer codes presented in Sect. 3.

3 Radiative-transfer codes

In the following subsections, each group provides a brief description of their code, highlighting the specific setup used in the calculations for this paper. Each code subsection follows a similar structure: brief description of the code (and whether it assumes a homologous velocity law); treatment of γ-ray energy deposition and non-thermal effects; computation of the temperature structure; treatment of excitation and ionisation; evaluation of the radiation field; sources of opacity and atomic data; resolution (spatial and frequency) and typical runtime. Table 2 summarises the physical ingredients and approximations used in each code for the test models considered here.

3.1 ARTIS

ARTIS² is a Monte Carlo radiative transfer code that uses the indivisible energy packet method of Lucy (2002). The code was originally described by Sim (2007) and Kromer & Sim (2009), with later improvements presented by Bulla et al. (2015) and Shingles et al. (2020). The code is three-dimensional and follows the time-evolution of the radiation field and state of the gas. It assumes a strictly homologous velocity law.

Injection of energy into the ejecta is calculated by following the deposition of γ-ray packets that are injected in accordance with the radioactive decays of ⁵⁶Ni and ⁵⁶Co, following Lucy (2005). Additional decay channels have been included in the studies of specific models. The simulated γ-ray transport is non-grey and takes into account Compton scattering, photoelectric absorption, and pair-creation opacities. In our standard runs, the code does not include the effects of excitation or ionisation by non-thermal particles. However, Shingles et al. (2020) presented updates to the code that include a Spencer-Fano treatment of non-thermal ionisation as required for late-phase modelling. Results obtained with this improved version (artisnebular) are included for late phases for the models in this study.

The electron temperature in each grid zone is estimated by balancing of heating and cooling rates (accounting for γ-ray and positron deposition, bound-bound, bound-free, and free-free processes). In its regular mode of operation (artis), the code uses an approximate non-LTE treatment to estimate the ionisation state in the ejecta (based on Monte Carlo photoionisation estimators; see Kromer & Sim 2009) and an LTE treatment of excitation is adopted. This approach has been used in most ofour published studies, and is also used in most of the artis calculations presented here. However, this method has limitations that become increasingly important at later phases (e.g. it neglects non-thermal heating and ionisation and tends to overestimate the plasma temperature at low densities due to incomplete treatment of cooling by forbidden lines). To improve these issues, Shingles et al. (2020) presented updates to the code that include a full non-LTE population solver (together with the Spencer-Fano solver mentioned above). Results obtained with this improved version (artisnebular) are included for late-phase calculations here.

Monte Carlo estimators are used to track the radiation field in each grid cell. In general, we use volume-based estimators (see Lucy 1999 or Noebauer & Sim 2019) to extract radiation-field-dependent quantities from the flight histories of our Monte Carlo quanta. In its standard mode of operation, artis uses detailed Monte Carlo estimators to obtain photoionisation rates from the radiation field, but relies on an estimated scaling for excited-state photoionisation and on a dilute black-body radiation field model for bound-bound excitation (see Kromer & Sim 2009 for details). However, the improved artisnebular version uses a more detailed frequency binned representation of the radiation field (see Shingles et al. 2020 for details). The code has the capacity to iterate on each time step with the aim of achieving consistency between the radiation field estimates and the packet transport in each step. However, in practise we find that this iteration is not required and we therefore generally simply use the radiation field quantities extracted from the previous time step to estimate the radiative rates that are needed for the current step.

In simulating ultraviolet to infrared photon transport, the code accounts for electron scattering, bound-bound, bound-free, and free-free processes. Bound-free and bound-bound processes are treated using the Macro Atom approach of Lucy (2002, 2003) and adopting the Sobolev approximation for line opacity. The code does not use an expansion opacity (or similar) but treats line opacity based on a frequency-ordered list of transitions treated in the Sobolev limit (i.e. no line overlap is taken into account).

In our simulations, atomic data are primarily drawn from the Kurucz atomic line lists (see Kromer & Sim 2009) – in Appendix A.1 we summarise the ions and numbers of levels and lines that we include. The photon transport is carried out on a 3D Cartesian grid that co-expands with the ejecta. The artis simulations included here were carried out on a 100³ grid. The resolution therefore corresponds to around 500–1000 km s⁻¹, depending on the model. The simulations are typically run on 1000 computer cores for 1–2 days.

3.2 CMFGEN

CMFGEN is a 1D, non-LTE, time-dependent radiative-transfer code that solves the transfer equation in the co-moving frame of spherical outflows. Details about the code, techniques, and approximations can be found in Hillier & Miller (1998), Hillier (2003, 2012), and (for SN calculations) in Hillier & Dessart (2012)³. The velocity law for the outflow is in general monotonic (but see e.g. Dessart et al. 2015 in the case of interacting SNe) and is assumed here to be homologous (such that ∂v/∂r = v/r).

In the present calculations, non-local energy deposition from radioactive decay is treated using a Monte-Carlo approach for γ-ray transport (Hillier & Dessart 2012). Non-thermal processes associated with the high-energy electrons produced by Compton scattering and photoelectric absorption of these γ rays are accounted for through a solution of the Spencer-Fano equation (Spencer & Fano 1954; Li et al. 2012).

The temperature structure is constrained through the energy equation that has the form:

$$\rho \frac{De}{Dt}-\frac{P}{\rho }\frac{D\rho }{Dt}=4\pi {\displaystyle \int \text{d}v\left({\chi }_v{J}_v-{\eta }_v\right)+{\dot{\epsilon }}_{\text{decay}},}$$(9)

where $\frac{D}{Dt}$ is the Lagrangian derivative, e is the internal energy per unit mass, P is the gas pressure, and ${\dot{\epsilon }}_{\text{decay}}$ is the energy absorbed per second per unit volume (see Hillier & Dessart 2012 for further details). We verify the accuracy of the solution by examining a global energy constraint (equivalent to conservation of flux in a static atmosphere; see Hillier & Dessart 2012 for details), and an equation describing energy conservation as applied to the heating and cooling of free electrons. These two equations are related to the above Eq. (9) via the transfer equation, and the rate equations, respectively (e.g. Hillier 2003). Because of model assumptions (e.g. the use of super-levels) these two equations are not satisfied exactly, but the errors (typical at the 1% level or smaller) are too small to affect the validity of the models. Processes contributing to electron heating and cooling include bound-free ionisation and recombination-collisional ionisation and recombination, collisional excitation and de-excitation, free-free emission, Auger ionisation, charge exchange reactions (primarily with H I and He I, and hence negligible in SN Ia ejecta), net cooling from non-thermal processes, and heating by radioactive decay.

Atom and ion-level populations are determined through a solution of the time-dependent rate equations, coupled to the above energy equation and the zeroth and first moments of the radiative-transfer equation (see below). We consider the following processes: bound-bound processes, bound-free processes-collisional ionisation and recombination, collisional excitation and de-excitation, Auger ionisation (Hillier 1987; Hillier & Miller 1998)’ and non-thermal excitation and ionisation (Li et al. 2012). We additionally consider further processes involving H and He (two-photon decay, charge-exchange reactions, and Rayleigh scattering), although these are negligible here (and completely absent from the toy models, which contain no H or He). To ease the solution of the rate equations, atomic levels are grouped into super levels (see Hillier & Miller 1998 for details).

The frequency-dependent mean intensity J_v is obtained via a solution of the time-dependent transfer equation in the co-moving frame to first order in v/c. More specifically, we solve the zeroth and first moment equations, which are closed using so-called Eddington factors f_v = K_v/J_v, where K_v is the second moment of the specific intensity (related to the radiation pressure). The Eddington factors are obtained from a formal solution of the time-independent transfer equation.

We consider the following sources of opacity: electron scattering, bound-free (including photoionisation from excited states), bound-bound⁴, free-free, and Auger ionisation. As noted earlier, Rayleigh scattering and two-photon processes (for H and He only) are also part of the global opacity budget but are negligible here.

A description of the sources of atomic data can be found in the Appendix (Sect. A.2). The number of levels (both super-levels and full levels) and corresponding number of bound, bound transitions are given in Tables A.2-A.8. We ignore weak transitions with a gf value⁵ below some cutoff, typically set to 10^-4. For the toy models, the following ions were included: Si II–IV, S II–IV, Ca II–IV, Fe I–V, Co II–VII, and Ni II–V. For the delayed-detonation DDC models the following ions were included: C I–IV, O I–IV, Ne I–III, Na I, Mg II–III, Al II–III, Si II–IV, S II–IV, Ar I–III, Ca II–IV, Sc II–III, Ti II–III, Cr II–IV, Mn II–III, Fe I–VII, Co II–VII, and Ni II–VII (we also include the ground states of Cl IV, K III, and V I for the sole purpose of tracking changes in the abundances of radioactive isotopes). For all the aforementioned ions, we also consider ionisations to and recombinations from the ground state of the next ionisation stage (e.g. Ni VIII in the case of Ni). As time proceeds and the temperature in the spectrum-formation region drops, the highest ionisation stages have a negligible impact on the RT solution. When this occurs, smaller model atoms are used for these ions, or the ions are removed altogether from the atomic model set.

The toy models were remapped onto a coarser spatial grid with 100 depth points. No remapping was performed for the DDC models. Typical wall-clock runtimes are of the order of 24h per time step on a single computing node with 8-12 CPUs, thus taking 2-3 months to complete for a sequence covering the first 200 days or so post explosion.

3.3 CRAB

CRAB is a 1D, implicit, Lagrangian radiation hydrodynamics code developed to study the light curves during SN outbursts and the corresponding outflows evolved from hydrostatic state up to homologous expansion (Utrobin 2004). Non-local energy deposition of γ rays from radioactive decay is determined by solving the corresponding one-group γ-ray transport with the approximation of an effective absorption opacity of 0.06 Y_e cm² g⁻¹. Positrons are assumed to deposit their energy locally. This energy deposition produces non-thermal ionisation and heating, the rates of which are taken from Kozma & Fransson (1992).

The radiation hydrodynamic equations include a time-dependent energy equation, which is based on the first law of thermodynamics and determines the gas temperature structure. In the outer, transparent and semitransparent layers of the SN ejecta, the local energy balance is in control of a net balance between heating and cooling processes, while in the inner, optically thick layers, it is determined by the diffusion of equilibrium radiation.

The code has two options for the treatment of atom and ion level populations. In option A, the elements H, He, C, N, O, Ne, Na, Mg, Si, S, Ar, Ca, and Fe, and the negative hydrogen ion H⁻ are included in the non-LTE ionisation balance. All elements but H are treated with three ionisation stages. The ionisation balance is controlled by the following elementary processes: photoionisation and radiative recombination, electron ionisation and three-particle recombination, and non-thermal ionisation. In option B, we use an LTE treatment of ionisation and excitation for elements from the option A list or all elements from H to Zn modelled with an arbitrary number of ionisation stages. Atomic and ionic level populations are determined using the Boltzmann formulae and the Saha equations for an element mixture with the local electron or non-equilibrium radiation temperature.

The time-dependent radiation transfer is treated in a one-group (grey) approximation in the outer, transparent and semitransparent layers of the SN envelope, while in the inner, optically thick layers where thermalisation of radiation takes place and LTE applies, the diffusion of equilibrium radiation is described in the approximation of radiative heat conduction. The bolometric luminosity of the SN is calculated by including retardation and limb-darkening effects.

In the inner, optically thick layers, the Rosseland mean opacity is evaluated for the local electron temperature, while in the outer, transparent and semitransparent layers non-LTE effects are taken into account when determining the mean opacities and the thermal emission coefficient. The mean opacities include processes of photoionisation, free-free absorption, Thomson scattering on free electrons, and Rayleigh scattering on neutral hydrogen. The contribution of lines to the expansion opacity is evaluated by the generalised formula of Castor et al. (1975) or by the formula of Blinnikov (1996) using the Sobolev approximation for line opacity. The expansion opacity in an expanding medium is treated with a thermalisation parameter set to 0.9 as recommended by Kozyreva et al. (2020) to model SNe Ia.

The sources of atomic data for processes in continuum can be found in (Utrobin 2004). Oscillator strengths of lines are taken from the Kurucz line database⁶ containing nearly 530000 lines. Energy level data are from the atomic spectra database of the National Institute of Standards and Technology.

The zoning of model toy06 with 808 zones is adequate for modelling a light curve. For this model, a typical runtime is of the order of 5 h on one CPU for the first 140 days.

3.4 KEPLER

KEPLER is a one-dimensional, implicit Lagrangian hydrodynamics code with appropriate physics for the study of massive stellar evolution and SNe (Weaver et al. 1978; Woosley et al. 2002). The radiation transport is flux-limited diffusion using a single temperature for the matter and radiation.

An important difference between KEPLER and some of the other codes used here is that KEPLER does not assume a coasting configuration. At early times in particular, when the matter is very optically thick, energy input from ⁵⁶Ni and ⁵⁶Co decay will both heat the matter and accelerate it. The correction to the kinetic energy is small for the cases studied here, but the integral of the emitted light will be less than the integral of the decay energy that is deposited (minus adiabatic losses).

Our approach to γ-ray energy deposition is discussed in the Appendix of Ensman & Woosley (1988). Since the early 1990s, a value of κ_γ = 0.054 cm⁻² g⁻¹ has been used for the global opacity parameter in KEPLER’s leakage scheme to model SNe Ia and that is the value used here. A better value can be obtained by comparing with Monte Carlo calculations for a given class of model, such as SN IIP, SN Ib, and so on. The value 0.054 cm² g⁻¹ is larger than the standard local grey opacity (0.06 ± 0.01)Y_e cm² g⁻¹ (e.g. Swartz et al. 1995), where Y_e ≈ 0.5, as it is used to calculate the effective column depth from the outer edge of a spherical zone to the surface. For all points except the geometric centre of the explosion, the angle-averaged column depth would be greater than along this radial line. The averaging is therefore approximated by taking a larger opacity. In reality, this number would vary with the radial distribution of ⁵⁶Ni and would be smaller if all the ⁵⁶Ni were at the centre.

The temperature structure is computed by solving the hydrodynamics equations including the effects of expansion and acceleration with energy input by radioactive decay and transported by radiative diffusion. Because of the way KEPLER handles flux-limited transport using a single temperature for the radiation and matter, the pressure in the outermost zones can be overestimated. This can lead to a small (of order 10%) overestimate of the conversion of radiation to kinetic energy in those zones at late times. To alleviate this problem, the luminosity is taken to be the maximum of the value at the (electron-scattering) photosphere and the zone that includes 95% of the mass. The former dominated the light curve until after peak.

An important parameter of the calculation is the atomic opacity used to transport thermal radiation. In KEPLER, the total ‘optical’ opacity consists of two parts: (a) electron scattering, which is calculated using a full Saha solve for all ionisation stages of all 19 species present; and (b) a constant additive opacity, κ_a, taken to reflect the effect of Doppler-broadened lines. The electron-scattering opacity is temperature-, density-, and composition dependent and therefore varies with location and time. The additive opacity is a constant everywhere for all time. Traditionally, we have used a value of 0.1 cm² g⁻¹ in our studies of Type Ia light curves, but the best value will depend on SN type. A comparison of SN Ib models calculated using CMFGEN and KEPLER (Ertl et al. 2020) suggested a value of 0.03 cm² g⁻¹ and we regard this as a lower limit for the average. The actual value varies with location and time in a complicated way. Here we adopt κ_a = 0.1 cm² g⁻¹ for the results presented in this paper.

KEPLER carries a nuclear network of 19 species (Weaver et al. 1978) which do not perfectly correspond to the species in the initial models provided. Care was taken to translate the ⁵⁶Ni and ⁵⁶Co abundances given to the zero-age ⁵⁶Ni and stable iron mass fractions which are used by KEPLER for energy generation. The species ¹²C and ¹⁶O were mapped without change. Other species such as ²⁰Ne, ²⁸Si, and so on were slightly augmented where necessary by adding in nearby odd-Z abundances; that is, ²⁰Ne included ²⁰Ne and ²³Na, ²⁴Mg included ²⁴Mg and ²⁷Al, and so on. As KEPLER does not compute spectra, this should have negligible consequences.

The four models were mapped into KEPLER. The total mass, kinetic energy, and ⁵⁶Ni mass were preserved to 0.05%. Because KEPLER is not a special relativistic code, zones with velocities greater than 0.1c were trimmed from the input. High-velocity zones would also cause difficulty with editing the luminosity at late times if the light-crossing time ceased to be negligible. Removing this high-velocity material decreased the kinetic energy of all models by about 2%. This should have negligible effects on the light curve. The zoning of the DDC models was relatively coarse by traditional KEPLER standards. After trimming the high-velocity zones, only 80 zones remained. The zoning of the toy models was better with roughly 600 zones remaining. Rezoning was not carried out. For a light curve with no nuclear burning, the zoning is adequate. All calculations took at most a few minutes on a laptop.

3.5 SEDONA

SEDONA is a time-dependent, multi-frequency Monte Carlo radiative transfer code originally developed to study SN light curves, spectra, and polarisation (Kasen et al. 2006). For this comparison study, the gas properties were tracked on 1D spherical Lagrangian zones and are evolved under the assumption of homologous expansion.

The radioactive decays of ⁵⁶Ni and ⁵⁶Co were tracked and used to source γ-ray packets. For the DDC10 model, radioactive decays of ⁴⁸Cr and⁴⁸ V were also included in addition to the ⁵⁶Ni decay chain. More detailed radioactive decay networks can be implemented in SEDONA, but only the previously mentioned radioactive-decay processes were used in this comparison study.

The γ-ray packets were transported using the Monte Carlo procedure subject to simplified treatments of bound-free absorption and Compton down-scattering. The γ-ray interactions were treated as heating terms that entered into the thermal heating balance that sourced a population of thermal photon packets.

The temperature in each shell was computed assuming radiative equilibrium. In more detail, the LTE assumption let us set the gas emissivity to the Planck function multiplied by the frequency-dependent line opacity. The temperature was then iteratively adjusted until the frequency-integrated emissivity was equal to the total radiative energy from both γ-rays and reemitted photon packets absorbed over the previous particle propagation step, as measured in the co-moving frame of the fluid.

For most of its published applications, SEDONA assumes LTE in order to compute opacities, although non-LTE capabilities have been developed and were recently applied (Roth et al. 2016; Shen et al. 2021). For this comparison study, only the LTE opacity mode was used in order to compute the excitation and ionisation states of the gas. More precisely, this means that ionisation fractions for each element were computed by simultaneously solving the Saha equation using a local temperature, and the charge conservation equation across all elements and isotopes present. Meanwhile, the bound-electron level populations within each ionisation stage were set by the Boltzmann factors given by the local temperature.

Photon packets were transported in three dimensions as a direct solution to the time-dependent radiative transfer equation. For all interactions with the gas, the packets were mapped to the 1D Lagrangian zones.

Once the level populations had been computed, the bound-bound opacity was computed using the expansion opacity formalism, as described by Eastman & Pinto (1993) and Kasen et al. (2006). For the toy01 and toy06 models, we only included the bound-bound opacities and electron scattering opacity, using the Thomson cross section. For DDC10, simplified bound-free and free-free opacities were also included in addition to the previously mentioned opacities, although they did not have a noticeable effect during the early stages of the explosion, which is what we wish to compare using these calculations.

A thermalisation parameter є of 1.0 was used for the toy01 and toy06 models. This means that all photon packets (other than γ-ray packets) that were absorbed were immediately thermalised, so that each absorption event was followed by re-emission of a photon packet with a frequency sampled from the thermal emissivity. For the DDC10 model, є was set to 0.8, so in that case 20% of absorbed photon packets (not including the γ-ray packets) were coherently scattered instead of having their frequency re-sampled.

The atomic data used for the bound-bound transitions were taken from Kurucz CD 1. This is a larger set of atomic data than CD 23. The details are described in Appendix A.3.

For this comparison study, the co-moving frequency grid for the thermalised photons (i.e. not the γ-rays) used 17 664 bins, with equal logarithmic spacing, ranging from 10¹⁴ Hz to 2 × 10¹⁶ Hz. The output spectra used 1044 frequency bins over the range 1.1 × 10¹⁴ Hz to 2 × 10¹⁶ Hz (~150 Å to ~2.7 µm). The time steps for the homologous expansion and radioactive decay began at approximately a few hours at the start of the calculation, and grew to no longer than 1 day, with a maximum time-step growth rate of 10%. One million γ-ray packets were emitted at each time step. As these γ rays deposited their energy, their packets were discarded, while one thermally sourced photon packet was emitted for each discarded γ-ray packet. With all of the settings described above, a calculation run to 100 days post-explosion required about 10 CPU hours, and could be efficiently performed in parallel across several hundred processors, reducing the elapsed wall time to less than 1 h.

3.6 STELLA

The multi-group radiation hydrodynamics code STELLA (Blinnikov & Bartunov 1993; Blinnikov et al. 1998, 2006) is capable of computing the evolution of the radiation field coupled to the hydrodynamics, as well as the bolometric light curve, spectral energy distribution, and resulting broad-band magnitudes and colours. Therefore, STELLA does not require the condition of homologous expansion and is capable of treating shock propagation and any dynamical processes taking place in the SN ejecta.

Energy deposition from ⁵⁶Ni and ⁵⁶Co radioactive decay is treated in a one-group diffusion approximation with an absorption opacity of 0.05 Y_e cm² g⁻¹ according to Swartz et al. (1995). The energy of positrons is thermalised locally.

The ionisation and level populations of a limited set of species (H, He, C, N, O, Ne, Na, Mg, Al, Si, S, Ar, Ca, stable Fe, stable Co, and stable Ni) is treated in the LTE approximation. Radiation is not treated in equilibrium with the matter. The colour temperature is estimated as a black body temperature via the least-squares method.

The opacity includes photoionisation, free-free absorption, and electron scattering processes assuming LTE for the plasma and line interactions. Radioactive ⁵⁶Ni and ⁵⁶Co contribute to the stable Fe when the line opacity is calculated. The expansion opacity is calculated according to (Eastman & Pinto 1993). The thermalisation parameter for the line opacity treatment is set to 0.9 (Kozyreva et al. 2020). The line opacity is calculated using a data base of 153 441 spectral lines partly from Kurucz & Bell (1995), Verner & Yakovlev (1995), and Verner et al. (1996b).

The spectral energy distribution is computed in the wavelength range from 1Å to 50 000 Å. The frequency range is divided into 129 bins with equal logarithmic spacing, in which the radiative transfer equations are solved at every time step. For the toy06 model, a higher resolution simulation with 629 frequency bins was run, which is the version used in the following sections for this model. The overall opacity might be slightly underestimated in the simulations with 629 frequency bins compared to simulations involving 129 bins, because the expansion opacity is calculated for the lines in the given bin and is not extended to another bin even if the velocity gradient is very high. The final (pseudo-)bolometric light curve is obtained by integrating the spectra over frequency.

To avoid numerical artefacts, ejecta models closer in time to the explosion were used (all are included in the data repository): 1 h post explosion for both toy models, and ~29 s and ~63 s post explosion for models DDC10 and DDC25, respectively. The toy models were computed on a lower-resolution grid (with 202 zones), while the DDC models were computed at a higher resolution (with 399 zones). The typical runtime on a single processor is a few hours at most.

3.7 SUMO

SUMO (Jerkstrand 2011; Jerkstrand et al. 2011, 2012) is a homologous non-LTE code with radiative transfer. It is specialised to calculate spectra and light curves in the post-peak phases of the SN. The code is written in Fortran and is parallelised with MPI.

Gamma-ray transfer is done by ray tracing using a grey opacity of 0.06Y_e cm² g⁻¹. Positrons can be transferred (using an effective opacity of 8.5 times the γ-ray one) but here they were assumed to be locally trapped. The cascade of non-thermal electrons following the scattering of gamma rays and positrons is solved for with the Spencer-Fano method (Kozma & Fransson 1992).

Zone temperatures are solved from the first law of thermodynamics considering heating by non-thermal processes and photoionisation, and cooling by net collisional line excitation, free-free emission, and recombination. The temperature is solved either in steady state (heating = cooling; Jerkstrand et al. 2012) or time-dependently (Pognan et al. 2022). Here, the steady-state mode was used.

The rate equations are solved by considering spontaneous radiative decay with Sobolev escape probabilities (assuming homology) modified to include continuum processes and strong line overlaps, thermal and non-thermal collisions (excitations and ionisations), photoexcitation and de-excitation, photoionisation, and recombination. Options exist to also include charge transfer processes and molecular chemistry (Liljegren et al. 2020), but these were not considered here. The rate equations can be solved either in steady state (inflow = outflow; Jerkstrand et al. 2012) or time dependently (Pognan et al. 2022). Here, the steady-state mode was used.

Radiative transfer is calculated with a Monte Carlo simulation that is iterated with the solvers for temperature, ionisation, and excitation. During the transfer, photon packets can experience electron scattering, bound-free and free-free absorption, and line absorption. The line transfer is resolved line by line (with Sobolev formalism) rather than using an expansion opacity formalism. Photoexcitation is either fully coupled to the non-LTE solutions (for low- and mid-lying levels), or decoupled (for high-lying levels), instead giving a fluorescence cascade on the spot. The radiation field is computed in steady state.

The atomic data come from a variety of sources, mostly described in Jerkstrand et al. (2011, 2012). Model atoms have LS-states resolved, that is they do not use super levels. In models computed up to 2022 (including the ones here), ions up to and including doubly ionised were included, although higher ions are now being included (Pognan et al. 2022). For the toy models computed here, the model atoms are Fe I (496 levels), Fe II (578 levels), Fe III (600 levels), Ni I (136 levels), Ni II (500 levels), Ni III (8 levels), Co I (317 levels), Co II (108 levels), Co III (306 levels), Si I (494 levels), Si II (77 levels), Si III (2 levels), S I (123 levels), S II (5 levels), S III (6 levels), Ca I (198 levels), Ca II (69 levels), and Ca III (1 level)⁷.

For the runs here, the ejecta were resampled to ∆v = 500 km s⁻¹, and truncated at v = 30 000 km s⁻¹. The radiative transfer was computed on a wavelength grid from 400 Å to 25 000 Å, with a logarithmic resolution dλ/λ = 10⁻³. Typical wall-clock runtimes for a single epoch are a few hours on a typical 128 core setup.

3.8 SuperNu

SuperNu is a multi-group LTE radiative-transfer code that employs Implicit Monte Carlo (IMC) and Discrete Diffusion Monte Carlo (DDMC; Wollaeger et al. 2013; Wollaeger & van Rossum 2014). IMC solves the thermal radiative-transfer equations semi-implicitly by treating some absorption and emission as instantaneous effective scattering (see e.g. Fleck & Cummings 1971). DDMC optimises IMC in optically thick regions of space (Densmore et al. 2012) and ranges of wavelength (Abdikamalov et al. 2012; Densmore et al. 2012) by replacing many low mean-free-path scattering events with single leakage events. SuperNu can apply IMC and DDMC in both static and homologous, semi-relativistically expanding atmospheres. The code has been verified by analytic and semi-analytic radiative transfer tests (Wollaeger et al. 2013) and on the W7 model of SNe Ia (Nomoto et al. 1984; Wollaeger & van Rossum 2014).

For the γ-ray transfer, SuperNu employs a constant absorption opacity of 0.06 Y_e cm² g⁻¹, where Y_e is ejecta gas electron fraction, following the prescription of Swartz et al. (1995). The γ-ray packets in SuperNu are not directly converted to optical packets, but instead are used to tally the total γ-ray energy deposition per spatial cell. The deposition energy values are then added to the thermal source for optical packets.

The ejecta gas temperature is calculated using the standard IMC semi-implicit linearisation (Fleck & Cummings 1971) of the comoving internal energy equation (Wollaeger et al. 2013): internal energy is recast to gas temperature using the standard relation ∂e = c_v∂T, where e is internal energy, c_v is heat capacity at constant volume, and T is gas temperature. Due to the IMC time linearisation, energy deposited from gamma-rays and (locally) from beta particles appears simultaneously in both the comoving internal energy equation and the radiation equation.

Ionisation and excitation are both treated with Saha–Boltzmann statistics evaluated at the gas density and temperature at the beginning of the time step. The multi-element system is solved iteratively by converging the free electron number. The resulting population densities are then used to calculate opacities. The radiation field is represented by fully time-dependent Monte Carlo packets, which are sourced from the LTE emissivity. The radiation field is not constrained to be in equilibrium with the gas, and so in general the system is ‘two-temperature’.

Opacity is discretised into groups via direct integration over co-moving wavelength (but see Fontes et al. 2020 for a study with SuperNu using expansion opacity). Opacity in SuperNu includes free-free (Sutherland 1998), and bound-free (Verner & Yakovlev 1995; Verner et al. 1996a) processes, as well as the bound-bound opacities from the Kurucz line lists⁸, and the standard elastic Thomson scattering opacity. Weak lines in the Kurucz data set are omitted from the SuperNu line list where the opacity is dominated by stronger lines. The total number of available lines for the present simulations is 591 288. This list was motivated by studies using the PHOENIX code in the work of van Rossum (2012), using the full line list (>10⁷ lines) as a benchmark. For the simulations in this work, opacity is computed in 1000 logarithmic wavelength groups from 100 to 32 000 Å.

Each simulation presented for SuperNu uses 4 194 304 source packets per time step, with maximum active packet populations of between 100 and 200 million. The wavelength group structure that the packets are tracked through is the same as the opacity group structure (though for the simulations here, the wavelength bounds for the flux tallies are 1000 and 32 000 Å). The DDC and toy models use 115 and 202 velocity space cells, respectively. For the time grid, the DDC and toy models use 200 logarithmically increasing time steps out to day 80, from about day 1 (DDC models) or day 2 (toy models) post explosion.

SuperNu has MPI+OpenMP parallelisation. The SuperNu simulations of the toy01, toy06, DDC10, and DDC25 models presented here cost 190, 320, 398, and 400 core-hours, respectively, each using 16 MPI ranks and 8 OpenMP threads per rank.

3.9 TARDIS

TARDIS is an open-source steady-state 1D radiative-transfer code that uses indivisible energy packets as its transport quanta following the methods in Abbott & Lucy (1985), Lucy & Abbott (1993), Mazzali & Lucy (1993), Lucy (2002, 2003, 2005). Kerzendorf & Sim (2014) describes the initial version of the code which was primarily used to model SNe Ia. Subsequently, the code has been significantly enhanced to include a non-thermal approximation treatment for helium (Boyle et al. 2017) and the continuum processes and relativity treatments required for hydrogen-rich SNe (Vogl et al. 2019). TARDIS has been continuously enhanced since then (see e.g. Kerzendorf et al. 2022, similarly gamma-ray energy deposition is a new module of TARDIS, but has not been used for the models in this work). Full documentation and an extended physics walk-through for TARDIS can be found online⁹.

TARDIS assumes a steady-state homologously expanding SN envelope and injects Monte Carlo packets – randomly sampled from a given distribution (by default a black-body) – from an inner boundary. The code supports bound-bound, bound-free, free-free, and Thomson opacities with several redistribution schemes from simple scattering to a macro-atom (Lucy 2002, 2003). Summary packet statistics are used to estimate radiation field quantities (temperature, dilution factors, mean intensities, heating- and photo-ionisation rates, and line source functions). The estimated quantities are used to calculate the ionisation and excitation populations in steady state with a choice between LTE and several formulations of non-LTE (nebular approximation of Abbott & Lucy 1985 for ionisation, and dilute-LTE and full non-LTE for excitation). TARDIS then calculates Sobolev optical depths for line interaction, and opacities for continuum processes. These values are then used in the subsequent Monte Carlo step, which produces summary statistics for updating the opacities. The other convergence criterion is a match of the integrated packet output luminosity and the requested luminosity, which is achieved through iterative adjustments of the inner temperature.

To handle atomic data, the TARDIS collaboration has developed an additional package named CARSUS (Pássaro et al. 2019)¹⁰ with documentation available online¹¹. CARSUS can read atomic data (masses, ionisation energies, levels, lines, photoionisation cross sections, and collisional cross sections) from NIST (Kramida & Ralchenko 1999), Chianti (Dere et al. 2019, 1997), CMFGEN (Hillier & Miller 1998; Dessart & Hillier 2010), and Kurucz (Kurucz 2009). For the models in this paper, we used Kurucz CD 23 as the source of atomic data (see Appendix A.4 for the full description).

The setup files for TARDIS that are used in this work are available online¹² using an atomic data set from Kurucz CD 23. In the setup used for the code comparison work, we run TARDIS in a mode that self-consistently finds a temperature stratification given an inner boundary velocity and output luminosity. We used the mean bolometric luminosity of the other comparison codes for our output luminosity. The structure of the toy01 and toy06 models was set at 520 shells between 9000 and 35 000 km s⁻¹. The DDC models were truncated to 40 shells between 9000 and 35 000 km s⁻¹. The inner boundary velocity is found by iterating until the dilution factor is close to 0.5 in the innermost zone. Table 3 shows the inner boundary velocity and requested output luminosity for the different epochs.

We used the nebular approximation for ionisation and the dilute-lte approximation for excitation. We used 5 × 10⁵ packets for estimating our radiation field and 30 iterations for convergence. In the final iteration, we estimated the source functions with 5 × 10⁵ packets and then used the formal integral to synthesise a spectrum. For the comparison, we use line opacities and Thomson opacities with the macroatom interaction scheme. The resolution of the output spectra was uniform from 500 to 20 000 Å with 2 Å bin width. The models were run on one CPU with runtimes of less than 1 h.

3.10 URILIGHT

URILIGHT is a time-dependent Monte-Carlo code written in Fortran 90 by Yoni Elbaz based on the approximations that are used in SEDONA (see Sect. 3.5 above, Kasen et al. 2006 and references therein), in particular assuming homologous expansion. A detailed description of this program and previous comparisons to other published radiative-transfer codes for several benchmark problems are presented in Wygoda et al. (2019)¹³.

Energy deposition resulting from the decay of radioactive isotopes is calculated by a Monte-Carlo solution of the γ-ray transport, for which interaction with matter is included through Compton scattering and photoelectric absorption. For the calculations in this work, only ⁵⁶Ni and ⁵⁶Co decay were included. The temperature is iteratively solved for in each cell by requiring that the total emissivity be equal to the total absorbed energy. LTE is assumed for calculating the ionisation and excitation states: ionisation is obtained by solving the Saha equation, and excitation levels are set by the Boltzmann-distribution.

Opacities include bound-bound and free-free absorption and Thomson scattering off free electrons. The atomic line data for the bound-bound transitions, which constitutes the main and most important source of opacity in SNe Ia, are taken from the extended set of lines by Kurucz¹⁴. Following a bound-bound interaction, most photons (a fraction є, which is a global parameter of the simulation; see Kasen et al. 2006) are thermalised and re-emitted at a different wavelength with a distribution set by the emissivity. The rest of the photons (fraction 1 − є) are re-emitted at the same wavelength (within the line width). As in Wygoda et al. (2019), we used є = 0.8 in the runs performed here.

The runs here use 162 spatial cells for the toy01 and toy06 models and 115 spatial cells for the DDC10 and DDC25 models. All models are run with 128 time steps logarithmically spaced between 2 and 210 days, and a uniform spectrum resolution of 10 Å between 100 and 30 000 Å. Each run typically took 30 h on a single core, or approximately 1 h when parallelised.

4 Data repository of test models and standardised outputs

The ejecta models and output files from the RT simulations of the different codes are provided in a new data repository, which is publicly available and can be accessed at online¹⁵.

Descriptions of the files available in the repository are provided below, including ejecta model files (Sect. 4.1), output files (Sect. 4.2), and Python codes that were used to create the analytic toy model ejecta and codes for reading the output files (Sect. 4.3).

4.1 Ejecta model files

As described in Sect. 2, the code comparison is performed using four SN Ia models (main parameters in Table 1). The RT input files that were distributed among the groups are provided in the repository, including two toy-model files snia_toy01_2d.dat and snia_toy06_2d.dat and the two delayed-detonation model files DDC10_0.976d.dat and DDC25_1.300d.dat.

4.1.1 Toy model files

The two toy-model files, snia_toy01_2d.dat and snia_toy06_2d.dat, represent a snapshot of the ejecta at 2.0 days post explosion and include 807 shells (rows) with the following 21 columns¹⁶:

(1) Shell index (1–807)

(2) Velocity at shell centre (km s⁻¹)

(3) Shell mass (M_⊙)

(4) Lagrangian mass coordinate at the outer shell boundary (M_⊙)

(5) Pre-decayed (t = 0) stable IGE mass fraction (0 for all shells)

(6) Pre-decayed (t = 0) ⁵⁶Ni mass fraction

(7) IME mass fraction (of which 10% is Ca, 35% is S, and 55% is Si by mass)

(8) Ti mass fraction (0 for all shells)

(9) Unburnt C+O mass fraction (0 for all shells)

(10) Radius at shell centre (cm) = velocity at shell centre × 2 days (homologous expansion)

(11) Mean density over shell (g cm⁻³), not density at shell centre

(12) Temperature (K)

(13)-(21) Mass fractions of ⁵⁶Ni, ⁵⁶Co, ⁵⁶Fe, Ca, S, Si, O, and C at 2 days post explosion

4.1.2 Delayed-detonation model files

The two delayed-detonation model files DDC10_0.976d.dat and DDC25_1.300d.dat represent a snapshot of the ejecta at 0.976 days and 1.3 days, respectively, with 115 shells with the following 50 columns¹⁷:

(1) Velocity at shell centre (km s⁻¹)

(2) Radius at shell centre (cm) = velocity at shell centre × the time post explosion

(3) Shell volume (cm³)

(4) Density at zone centre (g cm⁻³), not mean density over shell

(5) Shell mass estimate (g) = shell volume × density at zone centre

(6) Temperature (K)

(7)–(26) Elemental mass fractions at snapshot time of C, O, Ne, Na, Mg, Al, Si, S, Cl, Ar, K, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, and Ni

(27)–(50) Isotopic mass fractions of the following radioactive nuclei at snapshot time: ⁵⁶Ni, ⁵⁶Co, ⁵⁷Ni, ⁵⁷Co, ⁴⁸Cr, ⁴⁸V, ⁴⁹Cr, ⁴⁹V, ⁵¹Mn, ⁵¹Cr, ⁵⁵Co, ⁵⁵Fe, ³⁷K, ³⁷Ar, ⁵²Fe, ⁵²Mn, ⁴⁴Ti, ⁴⁴Sc,⁴¹ Ar,⁴² Ar, ⁴²K, ⁴³K, ⁴⁷Sc, and ⁶¹Co.

4.2 RT output files

For each simulation by a specific group, six output file types are generated for each of the ejecta models.

4.2.1 Output file names

The name of each file has the following structure:

<output type>_<model>_<code name>.txt

where

• <output type> represents the type of output (described below) and can take one of six values: lbol_edep, edep, phys, ionfrac_<element> (where <element> is the name of a given element, e.g. ca for calcium), spectra, and wsynphot_mags,

• <model> represents one of the models and can take one of four values: toy06, toy01, ddc10, and ddc25,

• <code name> represents the code name with an optional additional descriptor (useful to distinguish between different code settings when applied to a given model) and can take one of 12 values: artis, artisnebular, cmfgen, crab, kepler, sedona, stella, stella_fr600 (for the STELLA runs for the toy06 models that use 629 frequency bins instead of the default 129), supernu, sumo, tardis, and urilight.

In principle, there are 12 code names × 4 ejecta models × 5 output files (excluding the ionfrac_<element> files) = 240 files and an additional 12 code names × (2 toy models × 6 elements + 2 DDC models × 20 elements) = 624 ionfrac_<element> files. This results in a total of 864 files, although in practice not all files are available for various reasons: (a) a code was not applied to a given model (it was agreed that all groups should at least compute the toy06 model, but the other three were considered optional); (b) a given code cannot produce the specified output (e.g. SEDs for grey codes); (c) a given code does not provide the desired output quantities by default (i.e. modification of the source code would be necessary). Table 4 summarises the outputs and computed models for each code.

4.2.2 Six output file types

The six types of output files include:

1. Pseudobolometric (UVOIR) luminosity and global energy deposition as a function of time.

File name: lbol_edep_<model>_<code name>.txt

Header: after some optional comment lines, the following two lines give the number of epochs (NTIMES, 100 in this example) and the column headings:

#<optional comment lines>
#NTIMES: 100
#time[d] Lbol[erg/s] Edep[erg/s]

Contents: NTIMES rows with the following three columns:

(1) Time since explosion in days

(2) Pseudobolometric (UVOIR) luminosity in erg s⁻¹

(3) Global energy deposition by γ rays and positrons in erg s⁻¹.

2. Energy deposition

File name: edep_<model>_<code name>.txt

Header: after some optional comment lines, the following four lines give the number of epochs (NTIMES, 100 in this example), the number of cells (NVEL, 200 in this example), the list of epochs in days with all NTIMES values (the ‘…’ should correspond to actual values), and finally the column headings (here the ‘… ’ can be used as is):

#<optional comment lines>
#NTIMES: 100
#NVEL: 200
#TIMES[d]: 2.0 3.0 … 100.0
#vel_mid[km/s] Edep_t0[erg/s/cm^3] Edep_t1[erg/s/cm^3] … Edep_tn[erg/s/cm^3]

Contents: NVEL rows with the following NTIMES+1 columns (101 in this example):

(1) Velocity at the centre of each cell in km s⁻¹

(2)–(101) Total γ-ray + positron energy deposition rate at each ofthe NTIMES epochs within the corresponding cell in erg s⁻¹ cm⁻³.

3. Physical conditions

File name: phys_<model>_<code name>.txt

Header: after some optional comment lines, the following three lines give the number of epochs (NTIMES, 100 in this example), the list of epochs with all NTIMES values (the ‘… ’ should correspond to actual values), and finally one empty comment line:

#<optional comment lines>
#NTIMES: 100
#TIMES[d]: 2.0 3.0 … 100.0
#

Contents: NTIMES blocks (one for each epoch), each containing a block header with three lines giving the epoch (in days, 2.0 days in this example), the number of cells saved for this epoch (NVEL, 200 in this example), and finally the column headings:

#TIME: 2.0
#NVEL: 200
#vel_mid[km/s] temp[K] rho[gcc] ne[/cm^3] natom[/cm^3]

Each block content consists of NVEL rows with the following five columns:

(1) Velocity at the centre of each cell in km s⁻¹

(2) Temperature in K

(3) Density in g cm⁻³

(4) Free electron density in cm⁻³

(5) Total atom density in cm⁻³.

4. Ionisation fraction (one file per element)

File name:

ionfrac_<element>_<model>_<code name>.txt

Header: after some optional comment lines, the following four lines give the number of epochs (NTIMES, 100 in this example), the number of ionisation stages (NSTAGES, 6 in this example, starting at neutral and up to NSTAGES – 1 times ionised), the list of epochs with all NTIMES values (the ‘…’ should correspond to actual values), and finally one empty comment line:

#<optional comment lines>
#NTIMES: 100
#NSTAGES: 6
#TIMES[d]: 2.0 3.0 … 100.0
#

Contents: NTIMES blocks (one for each epoch), each containing a block header with three lines giving the epoch (in days, 2.0 days in this example), the number of cells saved for this epoch (NVEL, 200 in this example), and finally the column headings (where we consider the element Fe in this example):

#TIME: 2.0
#NVEL: 200
#vel_mid[km/s] fe0 fe1 fe2 fe3 fe4 fe5

Each block content consists of NVEL rows with the following NSTAGES+1 columns (7 in this example):

(1) Velocity at the centre of each cell in kms⁻¹

(2)–(7) Fraction of ions (dimensionless) in the corresponding cell (fe0 = Fe I, fe1 = Fe II etc.). The sum of the fractions in each of the NVEL rows is expected to be unity.

5. Spectral sequence

File name: spectra_<model>_<code name>.txt

Header: after some optional comment lines, the following four lines give the number of epochs (NTIMES, 100 in this example), the number of wavelengths (NWAVE, 2000 in this example), the list of epochs in days with all NTIMES values (the should correspond to actual values), and finally the column headings (here the ‘…’ can be used as is):

#<optional comment lines>
#NTIMES: 180
#NWAVE: 2000
#TIMES[d]: 2.0 3.0 … 100.0
#wavelength[Ang] flux_t0[erg/s/Ang]
flux_t1[erg/s/Ang] … flux_tn[erg/s/Ang]

Contents: NWAVE rows with the following NTIMES+1 columns (101 in this example):

(1) Wavelength in Å

(2)–(101) Fluxes at each of the NTIMES epochs at the corresponding wavelength in erg s⁻¹ Å⁻¹.

6. Synthetic photometry

File name: wsynphot_mags_<model>_<code name>.txt

Header: after some optional comment lines, the following three lines give the number of epochs (NTIMES, 100 in this example), the number of photometric bands (NBANDS, 8 in this example), and finally the column headings (giving each band name, UBVRIJHK in our example):

#<optional comment lines>
#NTIMES: 100
#NBANDS: 8
#time[d] U B V R I J H K

Contents: NTIMES rows with the following NBANDS+1 columns (9 in our example):

(1) Time since explosion in days

(2)–(9) Absolute magnitudes in each band at the corresponding time.

4.3 Useful Python codes

All the figures presented in this paper are automatically generated and accessible in a Python notebook all_plots.ipynb which is available in the data repository. A separate Python notebook photometry.ipynb in the code-comparison1/directory is used to generate the synthetic photometry (wsynphot_mags files) from the spectra files on the fly. Three further useful Python codes can be used to:

1. Generate the toy models

Location:

data1/input_models/mk_snia_toy_model.py

Calling syntax:

python mk_snia_toy_model.py --highni
python mk_snia_toy_model.py --lowni

This code generates the toy models used in this paper, as explained in Sect. 2.1. The upper command generates the toy06 model (snia_toy06_2d.dat), while the lower line generates the toy01 model (snia_toy01_2d.dat). Comments inside the code provide information on how to create new similar test models with different parameters.

2. Read the input files

Location: data1/input_models/read_inputs.py

Calling syntax:

python read_inputs.py

This code reads the input files of the test models and creates the following files corresponding to Figs. 1 and 2 as well as Table 1: density_profile.pdf (Fig. 1), composition_profile_<model>.pdf (where <model> is one of toy06, toy01, ddc10, and ddc25; Fig. 2), and models_summary.tex (Table 1). The code also includes the functions read_snia_toy_model() and read_ddc_model() to read the files into Python variables.

By default, the code is expected to be executed from within the data1/input_models/ directory, although the path to the data1/ directory can be specified using the --path2data option.

3. Read the output files

Location: data1/read_outputs.py

Calling syntax:

python read_outputs.py file1.txt file2.txt
python read_outputs.py /path/to/file*.txt

This code reads the six output file types (see 4.2.2) and produces corresponding plots. The code also includes the functions: read_lbol_edep() , read_spectra() , read_edep() , read_phys() , read_ionfrac() , and read_mags() to read the output files into Python variables.

The code can be executed from any directory because the full path of each file can be specified (wildcards are also accepted). When uploading new output files to the data repository users are required to ensure they match the expected format exactly. This can be achieved using the --checkformat option (and optionally the --noplot option to disable the plotting functionality when checking a large number of files).

Fig. 3 Pseudobolometric (UVOIR) light curves for model toy06. The inset shows a zoom into the maximum-light epoch (the estimated time of maximum light is indicated with a ‘+’ sign).

Fig. 4 Peak pseudobolometric (UVOIR) luminosities vs. rise times (from explosion to peak) for all four test models. The data points corresponding to the CMFGEN, CRAB, and SuperNu calculations of the toy06 model are difficult to distinguish.

5 Example results

In this section, we provide example results that are extracted from the outputs of the different simulations. The purpose of this section is to illustrate the contents of the data repository and comparisons that can be made using it. No comparison to observations or in-depth investigations of the sources of differences in results are made here. We note that while some groups provided a few sets of results for different physical approximations that appear in the repository, only one set of results is shown for each code in this set of examples. One exception is ARTIS, for which both the regular (artis) and nebular (artisnebular) versions are included for each model in order to present results for this code from early to late times.

The results are provided in the following subsections. In Sect. 5.1, bolometric properties of the emission are shown for simulations of the toy06 model, including the pseudobolometric (UVOIR) light curve and the energy deposition rate as a function of time. Rise times and peak luminosities are provided for the toy06 and toy01 models. In Sect. 5.2, the resulting profiles of temperature and mean ionisation state of cobalt are shown at selected times. In Sect. 5.3, multi-band light curves and colour curves are shown for model toy06. B-band rise times are provided for all models where available. Example spectra at selected early and late times are shown in Sect. 5.4. Throughout the section, each figure provides a legend with the colours or symbols associated with each code. A uniform coding scheme is used, such that a given code can be systematically identified in all the figures.

Fig. 5 Total energy deposition rate from y rays and positrons for model toy06, normalised to the analytic function given by Eq. (10) (the normalisation allows differences to be seen more clearly).

Fig. 6 Left: instantaneous ratio of the luminosity and the total energy deposition rate for model toy06. Right: ratio of the time-weighted integrals of the luminosity and energy deposition rate for model toy06 (here we do not show the results for the SUMO code because the calculation starts at 75 days post explosion). Both ratios are expected to reach unity at late times.

Fig. 7 Ejecta (gas) temperature and ionisation profiles as a function of velocity for the toy06 model. Left column: Temperature at 15 days (top), 40 days (middle), and 200 days (bottom) post explosion. Right column: mean ionisation of cobalt at 15 days (top), 40 days (middle), and 200 days (bottom) post explosion. We restrict the abscissa range of the ionisation plots to v ≤ 12 000 km s−1 since the Co abundance drops to 0 at larger velocities.

Fig. 8 Multi-band (UBVRIJHK) light curves for the toy06 model. We note the different ordinate range for the U band.

5.1 Bolometric (UVOIR) evolution

In this subsection, bolometric properties are shown based on the lbol_edep output files. In Fig. 3, the UVOIR light curves for model toy06 are shown for the first 100 days after explosion for the different codes. The rise times and peak luminosities are shown in Fig. 4 for all four test models. The peak luminosities were obtained by a parabolic fit to light curve around maximum light.

In Fig. 5, the total energy deposition from γ rays and positrons based on the same files is shown as a function of time for the first 200 days after explosion for the toy06 model. In order to highlight the differences among codes, all results are normalised to the same analytic approximation for the deposition:

$$\begin{array}{lll}{Q}_{\text{dep,norm}}\hfill & \equiv \hfill & \frac{M_{\text{Ni}}}{M_{\odot }}\left(0.97\left[1-e^{-{\left(40d/t\right)}^2}\right]+0.03\right)\hfill \\ \hfill & \hfill & \times \left(6.45e^{-t/8.8d}+1.45e^{-t/111.3d}\right){10}^{43}\text{erg\hspace{0.17em}}{\text{s}}^{-1}\hfill \end{array}$$(10)

where M_Ni = 0.6 M_⊙ is the total (undecayed) ⁵⁶Ni mass in the toy06 model¹⁸.

In Fig. 6, two physical diagnostics of the relation between the energy deposition and the bolometric light curves are shown for the toy06 model (both expected to approach unity at late times): the instantaneous ratio of the energy deposition rate and the luminosity (left) and the ratio of cumulative time-weighted integrals of the energy deposition rate and luminosity (right).

Fig. 9 Peak magnitude vs. rise time in the B band for all four test models.

Fig. 10 B – V and V – R colour evolution for the toy06 model.

Fig. 11 Spectra at 5 days post explosion and at the time of peak UVOIR luminosity (the caption indicates the time of the spectrum computed closest to peak) for the toy06 model.

5.2 Temperature and ionisation profiles

In Fig. 7, the thermodynamic structure of the ejecta as a function of velocity is shown based on the phys and ionfrac_co output files for model toy06. Three different times of particular interest are shown: 15 days post explosion, close to the peak of the light curve (upper panels), 40 days post explosion, close to the break in the B – V colour curve (middle panels), and 200 days post explosion, during the nebular phase (lower panels). For each of these times, two profiles are shown: the (gas) temperature profile (left panels), and the mean ionisation level of Co, defined as ∑_i=oⁱ fu where f_i is the fraction of Co ions ionised i times (right panels). Co and Fe dominate the opacity at these times and the ionisation profiles of Fe are very similar to those of Co (not shown here).

5.3 Multi-band light curves and colours

In this section, multi-band properties are shown for the different codes for different models. The photometry is extracted from the spectra reported in the spectra output files using the wsynphot package¹⁹ (using the Vega calibration spectrum alpha_lyr_stis_003. fits²⁰ and a third-order spline interpolation). In Fig. 8, the UBVRIJHK light curves are shown for model toy06 up to 100 days after explosion. The ß-band rise times (from explosion to peak) are extracted by fitting a high-order polynomial around maximum light and shown for all models in Fig. 9. The B – V and V – R colour curves for model toy06 are shown in Fig. 10.

5.4 Spectroscopic evolution

Spectra obtained for the toy06 model at selected times are shown in Figs. 11 and 12 based on the spectra output files. The wavelength range is restricted to 3000–10 000 Å to allow direct comparison. This does not represent the full range available in the output files, which differs between different codes. In Fig. 11, early-time spectra (5 days post explosion and around peak luminosity) are shown while in Fig. 12 we show spectra at later times (50 days and 200 days post explosion). We note that nebular times require specialised treatment in non-LTE which is only implemented in the ARTIS, CMFGEN, and SUMO codes for this paper²¹.

Fig. 12 Spectra at 50 days and 200 days post explosion for the toy06 model. The colour coding for the 50 days spectra is the same as in Fig. 11 and is not repeated here for sake of clarity.

6 Discussion and outlook

In this paper, we present an online public data repository that includes test models for radiative-transfer (RT) calculations of emission from SNe and standardised simulation outputs by ten different groups that allows direct comparison. Python scripts that generate the analytic toy models and read the different formats are provided as well.

The main purpose of the repository is to allow studies of the different physical approximations involved in the different codes and to assess the robustness of different predictions of radiation transfer. In addition, the repository can be used for finding bugs in the codes or can provide checks in the development of new codes.

We plan to extend the set of test models to include other SN types, as well as multi-dimensional models. In addition, we intend to produce more specialised test cases for which exact solutions can be found and agreed upon; these will provide benchmarks for RT code development. We also aim to include models with standardised atomic data sets that will allow the study of the effects of atomic physics on the emission for the different approximations.

Other groups that are developing RT codes are encouraged to add their results to the repository in the standardised format. Information for this purpose is provided in the README. md file in the repository.

Acknowledgements

The results presented in this paper are based on work performed before February 24, 2022. We thank the Schwartz/Reisman Institute for Theoretical Physics at the Weizmann institute of science and the Max Plank Institute for Astrophysics for hosting the workshops during 2018-2019 which lead to this collaboration. This work was supported by the ‘Programme National de Physique Stellaire’ (PNPS) of CNRS/INSU co-funded by CEA and CNES. This research was supported by the Excellence Cluster ORIGINS which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy EXC-2094-390783311. SB acknowledges support from the ESO Scientific Visitor Programme in Garching. CC acknowledges support by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 759253. A.J. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Program, ERC Starting Grant 803189 — SUPERSPEC. Lawrence Livermore National Laboratory is operated by Lawrence Livermore National Security, LLC, for the U.S. Department of Energy, National Nuclear Security Administration under Contract DE-AC52-07NA27344. AF and LS acknowledge support by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (ERC Advanced Grant KILO-NOVA No. 885281). The reported study was funded by the Russian Scientific Foundation (RSF), project number 19-12-00229, and the Russian Foundation for Basic Research (RFBR) and the Deutsche Forschungsgemeinschaft (DFG), project number 21-52-12032. S.A.S. acknowledges funding from the UKRI STFC Grant ST/T000198/1. Part of this work was performed using the Cambridge Service for Data Driven Discovery (CSD3), part of which is operated by the University of Cambridge Research Computing on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure. M.W. acknowledges support from the NASA Future Investigators in NASA Earth and Space Science and Technology grant (80NSSC21K1849) and support from the Thomas J. Moore Fellowship at New York University. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001). R.T.W. also acknowledges Daniel van Rossum for the development of SuperNu that enabled this work.

Appendix A Atomic data

In this Appendix we provide a more detailed overview of the atomic data used in the calculations done with the ARTIS (Sect. A.1), CMFGEN (Sect. A.2), SEDONA (Sect. A.3), and TARDIS (Sect. A.4) codes, along with appropriate references. We refer the reader to Sect. 3 for a more succinct presentation of the atomic data used by other groups and codes.

Appendix A.1 ARTIS

Table A.1 lists the elements with numbers of ions, levels, and bound-bound transitions used by the ARTIS code. These are drawn from the Kurucz atomic line lists (see Kromer & Sim 2009).

Appendix A.2 CMFGEN

Tables A.2-A.8 give the number of levels (both super-levels and full levels; see Hillier & Miller 1998 and Dessart & Hillier 2010 for details) for the model atoms used in the radiative-transfer calculations presented in this paper. N_SL refers to the number of super levels used for the solution of the rate equations, and NMI refer the number of full levels used to solve the transfer equation and compute the observed spectrum. We report the uppermost level for each ion treated in the fourth column. ‘W’ refers to states in which higher ℓ states (usually f or higher) have been combined into a single level. In the last column we give the number of bound-bound transitions in the model ion taking into account all N_full levels. Ions for which N_SL = N_full = 1 (with no entries for the last level configuration or number of lines) correspond to the final ionisation stage of a given element, for which ionisations to and recombinations from the ground state are considered. The ions Cl IV, KIII, and VI were included for the sole purpose of tracking changes in abundance of radioactive isotopes. The entries for those ions in Tables A.2-A.4 also have N_SL = N_fun = 1, but we specify the configuration of the ground state in the last level column and set the number of lines to zero.

Oscillator strengths for CNO elements were originally taken from Nussbaumer & Storey (1983, 1984). These authors also provide transition probabilities to states in the ion continuum. The largest source of oscillator data is from Kurucz (2009)²²; its principal advantage over many other sources (e.g. Opacity Project) is that LS coupling is not assumed. More recently, non-LS oscillator strengths have become available through the Iron Project (Hummer et al. 1993), and work done by the atomic-data group at Ohio State University (Nahar 2010). Other important sources of radiative data for Fe include Becker & Butler (1992, 1995a,b); Nahar (1995). Energy levels have generally been obtained from the National Institute of Standards and Technology. Collisional data are sparse, particularly for states far from the ground state. The principal source for collisional data among low-lying states for a variety of species is the tabulation by Mendoza (1983); other sources include Berrington et al. (1985), Lennon et al. (1985), Lennon & Burke (1994), Shine & Linsky (1974), Tayal (1997a,b), Zhang & Pradhan (1995a,b, 1997). Photoionisation data is taken from the Opacity Project (Seaton 1987) and the Iron Project (Hummer et al. 1993). Unfortunately Ni and Co photoionisation data are generally unavailable, and we used crude approximations. Charge exchange cross-sections are from the tabulation by Kingdon & Ferland (1996).

Appendix A.3 SEDONA

SEDONA used the Kurucz CD 1 line list to compute the bound-bound opacities.

Appendix A.4 TARDIS

TARDIS used the Kurucz CD 23 line list to compute the bound-bound opacities. The atomic data file was generated with the CARSUS package on August 24, 2017, named kurucz_cd23_chianti_H_He.h5 and signed with UUID1 6f7b09e887a311e7a06b246e9635S010 and MD5 864f1753714343c41f99cb06571Scace.

Table A.11 tabulates the total number of levels (N_levels), meta-stable levels (N_meta), and lines (N_lines) for the atoms used in the four models. Atoms present in the atomic data file but not used by the models were not listed.

References

All Tables

All Figures

	Fig. 1 Density profiles of our model set at an adopted time of 1 day post explosion.
In the text

Fig. 2 Composition profiles at the start of our simulations in our model set (2 days post explosion for the toy models and ~1 day post explosion for the DDC models). For the toy models, this represents the full set of species present (56Ni and decay products 56Co and 56Fe, as well as IMEs: Ca S, Si), while for the DDC models only a subset of species are shown for illustration. The 56Ni mass fraction is given at the time of explosion. We also show the total IGE mass fraction (from Sc to Ni; dashed line) and the total IME mass fraction (from Na to Ca; dotted line). The total IGE mass fraction coincides with the 56Ni0 line in the toy models and is not shown for sake of clarity.

In the text

	Fig. 3 Pseudobolometric (UVOIR) light curves for model toy06. The inset shows a zoom into the maximum-light epoch (the estimated time of maximum light is indicated with a ‘+’ sign).
In the text

	Fig. 4 Peak pseudobolometric (UVOIR) luminosities vs. rise times (from explosion to peak) for all four test models. The data points corresponding to the CMFGEN, CRAB, and SuperNu calculations of the toy06 model are difficult to distinguish.
In the text

	Fig. 5 Total energy deposition rate from y rays and positrons for model toy06, normalised to the analytic function given by Eq. (10) (the normalisation allows differences to be seen more clearly).
In the text

	Fig. 6 Left: instantaneous ratio of the luminosity and the total energy deposition rate for model toy06. Right: ratio of the time-weighted integrals of the luminosity and energy deposition rate for model toy06 (here we do not show the results for the SUMO code because the calculation starts at 75 days post explosion). Both ratios are expected to reach unity at late times.
In the text

Fig. 7 Ejecta (gas) temperature and ionisation profiles as a function of velocity for the toy06 model. Left column: Temperature at 15 days (top), 40 days (middle), and 200 days (bottom) post explosion. Right column: mean ionisation of cobalt at 15 days (top), 40 days (middle), and 200 days (bottom) post explosion. We restrict the abscissa range of the ionisation plots to v ≤ 12 000 km s−1 since the Co abundance drops to 0 at larger velocities.

In the text

	Fig. 8 Multi-band (UBVRIJHK) light curves for the toy06 model. We note the different ordinate range for the U band.
In the text

	Fig. 9 Peak magnitude vs. rise time in the B band for all four test models.
In the text

	Fig. 10 B – V and V – R colour evolution for the toy06 model.
In the text

	Fig. 11 Spectra at 5 days post explosion and at the time of peak UVOIR luminosity (the caption indicates the time of the spectrum computed closest to peak) for the toy06 model.
In the text

	Fig. 12 Spectra at 50 days and 200 days post explosion for the toy06 model. The colour coding for the 50 days spectra is the same as in Fig. 11 and is not repeated here for sake of clarity.
In the text