Issue 
A&A
Volume 668, December 2022



Article Number  A163  
Number of page(s)  25  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/202244134  
Published online  19 December 2022 
StaNdaRT: a repository of standardised test models and outputs for supernova radiative transfer
^{1}
Aix Marseille Univ, CNRS, CNES, LAM,
Marseille, France
email: stephane.blondin@lam.fr
^{2}
Unidad Mixta Internacional FrancoChilena de Astronomía, CNRS/INSU UMI 3386 and Instituto de Astrofísica, Pontificia Universidad Católica de Chile,
Santiago, Chile
^{3}
NRC “Kurchatov institute”, Institute for Theoretical and Experimental Physics (ITEP),
Moscow
117218, Russia
^{4}
Sternberg Astronomical Institute, Moscow State University,
Moscow
119234, Russia
^{5}
Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo,
515 Kashiwanoha, Kashiwa,
Chiba
2778583, Japan
^{6}
Astrophysics Research Centre, School of Mathematics and Physics, Queen’s University Belfast,
Belfast
BT7 1NN,
Northern Ireland, UK
^{7}
GSI Helmholtzzentrum für Schwerionenforschung,
Planckstraße 1,
64291
Darmstadt, Germany
^{8}
Institut d’Astrophysique de Paris, CNRSSorbonne Université,
98 bis boulevard Arago,
75014
Paris, France
^{9}
Center for Theoretical Astrophysics, Los Alamos National Laboratory,
Los Alamos, NM
87545, USA
^{10}
Theoretical Division, Los Alamos National Laboratory,
Los Alamos, NM,
87545, USA
^{11}
Department of Physics and Astronomy, Michigan State University,
East Lansing, MI
48823, USA
^{12}
Department of Physics and Astronomy & Pittsburgh Particle Physics, Astrophysics, and Cosmology Center (PITT PACC), University of Pittsburgh,
3941 O’Hara Street,
Pittsburgh, PA
15260, USA
^{13}
The Oskar Klein Centre, Department of Astronomy, Stockholm University,
AlbaNova,
SE10691
Stockholm, Sweden
^{14}
Departments of Physics and Astronomy, University of California Berkeley and Lawrence Berkeley National Laboratory,
USA
^{15}
Department of Particle Phys. & Astrophys., Weizmann Institute of Science,
Rehovot
76100, Israel
^{16}
Department of Computational Mathematics, Science, and Engineering, Michigan State University,
East Lansing, MI
48824, USA
^{17}
MaxPlanckInstitut für Astrophysik,
KarlSchwarzschildStraße 1,
85748
Garching bei München, Germany
^{18}
Facultad de Ciencias Astronómicas y Geofísicas, Universidad Nacional de La Plata,
La Plata
B1900, Argentina
^{19}
Department of Astronomy and Theoretical Astrophysics Center, University of California,
Berkeley, CA
94720, USA
^{20}
Institute of Astronomy, Russian Academy of Sciences,
Pyatnitskaya St. 48,
119017
Moscow, Russia
^{21}
Department of Physics, New York University,
New York, NY
10003, USA
^{22}
Computer, Computational, and Statistical Sciences Division, Los Alamos National Laboratory,
Los Alamos, NM
87545, USA
^{23}
Department of Physics & Astronomy, Louisiana State University,
Baton Rouge, LA
70803, USA
^{24}
Department of Astronomy and Astrophysics, University of California,
Santa Cruz, CA
95064, USA
^{25}
Department of Physics, NRCN,
BeerSheva
84190, Israel
^{26}
ARTIS Collaboration
^{27}
CMFGEN Collaboration
^{28}
CRAB Collaboration
^{29}
KEPLER Collaboration
^{30}
SEDONA Collaboration
^{31}
STELLA Collaboration
^{32}
SUMO Collaboration
^{33}
SuperNu Collaboration
^{34}
TARDIS Collaboration
^{35}
URILIGHT Collaboration
Received:
27
May
2022
Accepted:
22
September
2022
We present the first results of a comprehensive supernova (SN) radiativetransfer (RT) codecomparison initiative (StaNdaRT), where the emission from the same set of standardised test models is simulated by currently used RT codes. We ran a total of ten codes on a set of four benchmark ejecta models of Type Ia SNe. We consider two subChandrasekharmass (M_{tot} = 1.0 M_{⊙}) toy models with analytic density and composition profiles and two Chandrasekharmass delayeddetonation models that are outcomes of hydrodynamical simulations. We adopt spherical symmetry for all four models. The results of the different codes, including the light curves, spectra, and the evolution of several physical properties as a function of radius and time are provided in electronic form in a standard format via a public repository. We also include the detailed test model profiles and several Python scripts for accessing and presenting the input and output files. We also provide the code used to generate the toy models studied here. In this paper, we describe the test models, radiativetransfer codes, and output formats in detail, and provide access to the repository. We present example results of several key diagnostic features.
Key words: supernovae: general / radiative transfer
© S. Blondin et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
Accurate radiativetransfer (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 backreaction 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 codecomparison 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^{1}. 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 subChandrasekharmass (subM_{Ch}) toy models with profiles that are defined analytically and two more realistic models that result from hydrodynamical simulations of the M_{Ch} delayeddetonation 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.
Summary of ejecta conditions.
2 Test models
The codecomparison test suite consists of four SN Ia models. Two are subM_{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 ^{56}Ni yield, in order to have two models corresponding to ‘normal’ SNe Ia (toy06 and DDC10 with ~0.6 M_{⊙} of ^{56}Ni) and two lowluminosity models (toy01 and DDC25 with ~0.1 M_{⊙} of ^{56}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 lowluminosity SN Ia model) options. Both models have a total mass M_{ej} = 1.0 M_{⊙}, a kinetic energy E_{kin} = 10^{51} erg, and are calculated at the time t_{f} = 2 days post explosion. The models have an exponential density profile (e.g. Jeffery 1999),
where
is the efolding velocity, and
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^{−1}. 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:
where the integral runs from z_{0,i} to z_{1,i} and the volume element is , where we assume a homologously expanding ejecta (r = vt = v_{e}zt).
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:
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^{−1}. 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 irongroup elements (IGEs), here only ^{56}Ni and its decay products ^{56}Co and ^{56}Fe, and an outer region composed of the intermediatemass 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 delayeddetonation model DDC10 (see Sect. 2.2 below). Our ‘normal’ toy model consists of 0.6 M_{⊙} of ^{56}Ni and 0.4 M_{⊙} of IMEs, while our lowluminosity toy model consists of 0.1 M_{⊙} of ^{56}Ni and 0.9 M_{⊙} of IMEs. The ^{56}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 lowluminosity model and 0.4 M_{⊙} for the regular model). At a given mass coordinate m, the ^{56}Ni mass fraction is set to:
where m_{1} = M(^{56}Ni) − ∆M_{trans}/2, m_{2} = M(^{56}Ni) + ∆M_{trans}/2, and
The IME mass fraction is then simply set to . Our toy models therefore consist of only six chemical species or isotopes (^{56}Ni, ^{56}Co, ^{56}Fe, ^{40}Ca, ^{32}S, ^{28}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 ^{56}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 timeweighted internal energy, tE(t), is equal to the timeintegrated timeweighted decay energy deposition rate, ∫tQ(t)dt, as noted by Katz et al. (2013):
where X_{i}(^{56}Ni)_{0} is the ^{56}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 ^{56}Ni→^{56}Co→^{56}Fe decay chain. In this formulation, we ignore the timeweighted internal energy shortly after explosion, E(t_{0})t_{0}.
It is clear from Eq. (8) that the temperature is predicted to be zero in zones devoid of ^{56}Ni (≳12 000km s^{−1}). 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 longerterm 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 (^{56}Ni and decay products ^{56}Co and ^{56}Fe, as well as IMEs: Ca S, Si), while for the DDC models only a subset of species are shown for illustration. The ^{56}Ni 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 ^{56}Ni_{0} line in the toy models and is not shown for sake of clarity. 
2.2 Delayeddetonation models
In addition to the two toy models above, we consider two M_{Ch} delayeddetonation 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 ^{56}Ni yields of the toy models: Our lowluminosity model DDC25 yields ~0.12 M_{⊙} of ^{56}Ni (cf. 0.1 M_{⊙} for the toy01 model), and our ‘normal’ DDC10 model yields ~0.52 M_{⊙} of ^{56}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^{−1} 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 twostep ^{56}Ni→^{56}Co→^{56}Fe decay chain, we also treat eight additional twostep decay chains associated with ^{37}K, ^{44}Ti, ^{48}Cr, ^{49}Cr,^{51} Mn, ^{52}Fe, ^{55}Co, and ^{57}Ni, and a further six onestep decay chains associated with ^{41}Ar, ^{42}K, ^{43}K, ^{43}Sc, ^{47}Sc, and ^{61}Co (see Dessart et al. 2014).
The ejecta at 0.5 day are then remapped onto the 1D, nonLTE, radiativetransfer 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 radiativetransfer codes presented in Sect. 3.
3 Radiativetransfer 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 nonthermal 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^{2} 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 threedimensional and follows the timeevolution 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 ^{56}Ni and ^{56}Co, following Lucy (2005). Additional decay channels have been included in the studies of specific models. The simulated γray transport is nongrey and takes into account Compton scattering, photoelectric absorption, and paircreation opacities. In our standard runs, the code does not include the effects of excitation or ionisation by nonthermal particles. However, Shingles et al. (2020) presented updates to the code that include a SpencerFano treatment of nonthermal ionisation as required for latephase 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, boundbound, boundfree, and freefree processes). In its regular mode of operation (artis), the code uses an approximate nonLTE 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 nonthermal 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 nonLTE population solver (together with the SpencerFano solver mentioned above). Results obtained with this improved version (artisnebular) are included for latephase calculations here.
Monte Carlo estimators are used to track the radiation field in each grid cell. In general, we use volumebased estimators (see Lucy 1999 or Noebauer & Sim 2019) to extract radiationfielddependent 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 excitedstate photoionisation and on a dilute blackbody radiation field model for boundbound 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, boundbound, boundfree, and freefree processes. Boundfree and boundbound 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 frequencyordered 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 coexpands with the ejecta. The artis simulations included here were carried out on a 100^{3} grid. The resolution therefore corresponds to around 500–1000 km s^{−1}, depending on the model. The simulations are typically run on 1000 computer cores for 1–2 days.
Physical ingredients and approximations used in each code for the test models in this paper.
3.2 CMFGEN
CMFGEN is a 1D, nonLTE, timedependent radiativetransfer code that solves the transfer equation in the comoving 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)^{3}. 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, nonlocal energy deposition from radioactive decay is treated using a MonteCarlo approach for γray transport (Hillier & Dessart 2012). Nonthermal processes associated with the highenergy electrons produced by Compton scattering and photoelectric absorption of these γ rays are accounted for through a solution of the SpencerFano equation (Spencer & Fano 1954; Li et al. 2012).
The temperature structure is constrained through the energy equation that has the form:
where is the Lagrangian derivative, e is the internal energy per unit mass, P is the gas pressure, and 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 superlevels) 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 boundfree ionisation and recombinationcollisional ionisation and recombination, collisional excitation and deexcitation, freefree emission, Auger ionisation, charge exchange reactions (primarily with H I and He I, and hence negligible in SN Ia ejecta), net cooling from nonthermal processes, and heating by radioactive decay.
Atom and ionlevel populations are determined through a solution of the timedependent rate equations, coupled to the above energy equation and the zeroth and first moments of the radiativetransfer equation (see below). We consider the following processes: boundbound processes, boundfree processescollisional ionisation and recombination, collisional excitation and deexcitation, Auger ionisation (Hillier 1987; Hillier & Miller 1998)’ and nonthermal excitation and ionisation (Li et al. 2012). We additionally consider further processes involving H and He (twophoton decay, chargeexchange 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 frequencydependent mean intensity J_{v} is obtained via a solution of the timedependent transfer equation in the comoving frame to first order in v/c. More specifically, we solve the zeroth and first moment equations, which are closed using socalled 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 timeindependent transfer equation.
We consider the following sources of opacity: electron scattering, boundfree (including photoionisation from excited states), boundbound^{4}, freefree, and Auger ionisation. As noted earlier, Rayleigh scattering and twophoton 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 superlevels and full levels) and corresponding number of bound, bound transitions are given in Tables A.2A.8. We ignore weak transitions with a gf value^{5} 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 delayeddetonation 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 spectrumformation 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 wallclock runtimes are of the order of 24h per time step on a single computing node with 812 CPUs, thus taking 23 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). Nonlocal energy deposition of γ rays from radioactive decay is determined by solving the corresponding onegroup γray transport with the approximation of an effective absorption opacity of 0.06 Y_{e} cm^{2} g^{−1}. Positrons are assumed to deposit their energy locally. This energy deposition produces nonthermal ionisation and heating, the rates of which are taken from Kozma & Fransson (1992).
The radiation hydrodynamic equations include a timedependent 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 nonLTE 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 threeparticle recombination, and nonthermal 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 nonequilibrium radiation temperature.
The timedependent radiation transfer is treated in a onegroup (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 limbdarkening 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 nonLTE effects are taken into account when determining the mean opacities and the thermal emission coefficient. The mean opacities include processes of photoionisation, freefree 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^{6} 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 onedimensional, 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 fluxlimited 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 ^{56}Ni and ^{56}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^{−2} g^{−1} 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^{2} g^{−1} is larger than the standard local grey opacity (0.06 ± 0.01)Y_{e} cm^{2} g^{−1} (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 angleaveraged 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 ^{56}Ni and would be smaller if all the ^{56}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 fluxlimited 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 (electronscattering) 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 Dopplerbroadened lines. The electronscattering 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^{2} g^{−1} 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^{2} g^{−1} 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^{2} g^{−1} 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 ^{56}Ni and ^{56}Co abundances given to the zeroage ^{56}Ni and stable iron mass fractions which are used by KEPLER for energy generation. The species ^{12}C and ^{16}O were mapped without change. Other species such as ^{20}Ne, ^{28}Si, and so on were slightly augmented where necessary by adding in nearby oddZ abundances; that is, ^{20}Ne included ^{20}Ne and ^{23}Na, ^{24}Mg included ^{24}Mg and ^{27}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 ^{56}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. Highvelocity zones would also cause difficulty with editing the luminosity at late times if the lightcrossing time ceased to be negligible. Removing this highvelocity 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 highvelocity 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 timedependent, multifrequency 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 ^{56}Ni and ^{56}Co were tracked and used to source γray packets. For the DDC10 model, radioactive decays of ^{48}Cr and^{48} V were also included in addition to the ^{56}Ni decay chain. More detailed radioactive decay networks can be implemented in SEDONA, but only the previously mentioned radioactivedecay processes were used in this comparison study.
The γray packets were transported using the Monte Carlo procedure subject to simplified treatments of boundfree absorption and Compton downscattering. 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 frequencydependent line opacity. The temperature was then iteratively adjusted until the frequencyintegrated 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 comoving frame of the fluid.
For most of its published applications, SEDONA assumes LTE in order to compute opacities, although nonLTE 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 boundelectron 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 timedependent 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 boundbound 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 boundbound opacities and electron scattering opacity, using the Thomson cross section. For DDC10, simplified boundfree and freefree 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 reemission 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 resampled.
The atomic data used for the boundbound 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 comoving frequency grid for the thermalised photons (i.e. not the γrays) used 17 664 bins, with equal logarithmic spacing, ranging from 10^{14} Hz to 2 × 10^{16} Hz. The output spectra used 1044 frequency bins over the range 1.1 × 10^{14} Hz to 2 × 10^{16} 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 timestep 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 postexplosion 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 multigroup 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 broadband 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 ^{56}Ni and ^{56}Co radioactive decay is treated in a onegroup diffusion approximation with an absorption opacity of 0.05 Y_{e} cm^{2} g^{−1} 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 leastsquares method.
The opacity includes photoionisation, freefree absorption, and electron scattering processes assuming LTE for the plasma and line interactions. Radioactive ^{56}Ni and ^{56}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 lowerresolution 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 nonLTE code with radiative transfer. It is specialised to calculate spectra and light curves in the postpeak phases of the SN. The code is written in Fortran and is parallelised with MPI.
Gammaray transfer is done by ray tracing using a grey opacity of 0.06Y_{e} cm^{2} g^{−1}. 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 nonthermal electrons following the scattering of gamma rays and positrons is solved for with the SpencerFano method (Kozma & Fransson 1992).
Zone temperatures are solved from the first law of thermodynamics considering heating by nonthermal processes and photoionisation, and cooling by net collisional line excitation, freefree emission, and recombination. The temperature is solved either in steady state (heating = cooling; Jerkstrand et al. 2012) or timedependently (Pognan et al. 2022). Here, the steadystate 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 nonthermal collisions (excitations and ionisations), photoexcitation and deexcitation, 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 steadystate 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, boundfree and freefree 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 nonLTE solutions (for low and midlying levels), or decoupled (for highlying 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 LSstates 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)^{7}.
For the runs here, the ejecta were resampled to ∆v = 500 km s^{−1}, and truncated at v = 30 000 km s^{−1}. The radiative transfer was computed on a wavelength grid from 400 Å to 25 000 Å, with a logarithmic resolution dλ/λ = 10^{−3}. Typical wallclock runtimes for a single epoch are a few hours on a typical 128 core setup.
3.8 SuperNu
SuperNu is a multigroup LTE radiativetransfer 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 radiativetransfer equations semiimplicitly 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 meanfreepath scattering events with single leakage events. SuperNu can apply IMC and DDMC in both static and homologous, semirelativistically expanding atmospheres. The code has been verified by analytic and semianalytic 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^{2} g^{−1}, 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 semiimplicit 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 gammarays 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 multielement 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 timedependent 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 ‘twotemperature’.
Opacity is discretised into groups via direct integration over comoving wavelength (but see Fontes et al. 2020 for a study with SuperNu using expansion opacity). Opacity in SuperNu includes freefree (Sutherland 1998), and boundfree (Verner & Yakovlev 1995; Verner et al. 1996a) processes, as well as the boundbound opacities from the Kurucz line lists^{8}, 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^{7} 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 corehours, respectively, each using 16 MPI ranks and 8 OpenMP threads per rank.
3.9 TARDIS
TARDIS is an opensource steadystate 1D radiativetransfer 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 nonthermal approximation treatment for helium (Boyle et al. 2017) and the continuum processes and relativity treatments required for hydrogenrich SNe (Vogl et al. 2019). TARDIS has been continuously enhanced since then (see e.g. Kerzendorf et al. 2022, similarly gammaray 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 walkthrough for TARDIS can be found online^{9}.
TARDIS assumes a steadystate homologously expanding SN envelope and injects Monte Carlo packets – randomly sampled from a given distribution (by default a blackbody) – from an inner boundary. The code supports boundbound, boundfree, freefree, and Thomson opacities with several redistribution schemes from simple scattering to a macroatom (Lucy 2002, 2003). Summary packet statistics are used to estimate radiation field quantities (temperature, dilution factors, mean intensities, heating and photoionisation 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 nonLTE (nebular approximation of Abbott & Lucy 1985 for ionisation, and diluteLTE and full nonLTE 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)^{10} with documentation available online^{11}. 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^{12} 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 selfconsistently 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^{−1}. The DDC models were truncated to 40 shells between 9000 and 35 000 km s^{−1}. 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 dilutelte approximation for excitation. We used 5 × 10^{5} packets for estimating our radiation field and 30 iterations for convergence. In the final iteration, we estimated the source functions with 5 × 10^{5} 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.
Velocity and temperature at the inner boundary given the requested output luminosity for the TARDIS calculation of the toy06 model.
3.10 URILIGHT
URILIGHT is a timedependent MonteCarlo 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 radiativetransfer codes for several benchmark problems are presented in Wygoda et al. (2019)^{13}.
Energy deposition resulting from the decay of radioactive isotopes is calculated by a MonteCarlo 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 ^{56}Ni and ^{56}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 Boltzmanndistribution.
Opacities include boundbound and freefree absorption and Thomson scattering off free electrons. The atomic line data for the boundbound transitions, which constitutes the main and most important source of opacity in SNe Ia, are taken from the extended set of lines by Kurucz^{14}. Following a boundbound interaction, most photons (a fraction є, which is a global parameter of the simulation; see Kasen et al. 2006) are thermalised and reemitted at a different wavelength with a distribution set by the emissivity. The rest of the photons (fraction 1 − є) are reemitted 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^{15}.
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 toymodel files snia_toy01_2d.dat and snia_toy06_2d.dat and the two delayeddetonation model files DDC10_0.976d.dat and DDC25_1.300d.dat.
4.1.1 Toy model files
The two toymodel 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^{16}:
(1) Shell index (1–807)
(2) Velocity at shell centre (km s^{−1})
(3) Shell mass (M_{⊙})
(4) Lagrangian mass coordinate at the outer shell boundary (M_{⊙})
(5) Predecayed (t = 0) stable IGE mass fraction (0 for all shells)
(6) Predecayed (t = 0) ^{56}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^{−3}), not density at shell centre
(12) Temperature (K)
(13)(21) Mass fractions of ^{56}Ni, ^{56}Co, ^{56}Fe, Ca, S, Si, O, and C at 2 days post explosion
4.1.2 Delayeddetonation model files
The two delayeddetonation 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^{17}:
(1) Velocity at shell centre (km s^{−1})
(2) Radius at shell centre (cm) = velocity at shell centre × the time post explosion
(3) Shell volume (cm^{3})
(4) Density at zone centre (g cm^{−3}), 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: ^{56}Ni, ^{56}Co, ^{57}Ni, ^{57}Co, ^{48}Cr, ^{48}V, ^{49}Cr, ^{49}V, ^{51}Mn, ^{51}Cr, ^{55}Co, ^{55}Fe, ^{37}K, ^{37}Ar, ^{52}Fe, ^{52}Mn, ^{44}Ti, ^{44}Sc,^{41} Ar,^{42} Ar, ^{42}K, ^{43}K, ^{47}Sc, and ^{61}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.
Code outputs and computed models.
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^{−1}
(3) Global energy deposition by γ rays and positrons in erg s^{−1}.
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^{−1}
(2)–(101) Total γray + positron energy deposition rate at each ofthe NTIMES epochs within the corresponding cell in erg s^{−1} cm^{−3}.
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^{−1}
(2) Temperature in K
(3) Density in g cm^{−3}
(4) Free electron density in cm^{−3}
(5) Total atom density in cm^{−3}.
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^{−1}
(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^{−1} Å^{−1}.
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 codecomparison1/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 maximumlight 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 indepth 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, multiband light curves and colour curves are shown for model toy06. Bband 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 timeweighted 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 Multiband (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:
where M_{Ni} = 0.6 M_{⊙} is the total (undecayed) ^{56}Ni mass in the toy06 model^{18}.
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 timeweighted 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^{i} 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 Multiband light curves and colours
In this section, multiband 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^{19} (using the Vega calibration spectrum alpha_lyr_stis_003. fits^{20} and a thirdorder 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 highorder 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, earlytime 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 nonLTE which is only implemented in the ARTIS, CMFGEN, and SUMO codes for this paper^{21}.
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 radiativetransfer (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 multidimensional 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 20182019 which lead to this collaboration. This work was supported by the ‘Programme National de Physique Stellaire’ (PNPS) of CNRS/INSU cofunded 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 EXC2094390783311. 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 DEAC5207NA27344. 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 KILONOVA No. 885281). The reported study was funded by the Russian Scientific Foundation (RSF), project number 191200229, and the Russian Foundation for Basic Research (RFBR) and the Deutsche Forschungsgemeinschaft (DFG), project number 215212032. 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 eInfrastructure. 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 boundbound 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.2A.8 give the number of levels (both superlevels and full levels; see Hillier & Miller 1998 and Dessart & Hillier 2010 for details) for the model atoms used in the radiativetransfer 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 boundbound 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.2A.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)^{22}; its principal advantage over many other sources (e.g. Opacity Project) is that LS coupling is not assumed. More recently, nonLS oscillator strengths have become available through the Iron Project (Hummer et al. 1993), and work done by the atomicdata 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 lowlying 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 crosssections are from the tabulation by Kingdon & Ferland (1996).
Summary of atomic data used in ARTIS simulations.
Model atoms used in CMFGEN calculations for the toy06 model.
Model atoms used in CMFGEN calculations for the toy01 model.
Model atoms used in CMFGEN calculations up until 10.56 d post explosion for the DDC10 and DDC25 models.
Large Co model atoms used in CMFGEN calculations for the DDC10 (between 11.62d and 15.47 d post explosion) and DDC25 (between 11.62dand 17.02 d post explosion) models.
Model atoms used in CMFGEN calculations from 17.02 d post explosion for the DDC10 model and from 18.72 d post explosion for the DDC25 model.
Appendix A.3 SEDONA
SEDONA used the Kurucz CD 1 line list to compute the boundbound opacities.
Large Co model atoms used in CMFGEN calculations for the DDC10 and DDC25 models from 33.15 d post explosion onward.
Reduced Fe IV, Co IV, and Ni IV model atoms used in CMFGEN calculations for the DDC10 and DDC25 models from 40.11 d post explosion onward.
Boundbound atomic data information from the Kurucz CD 1 line list for the toy06 and toy01 models computed with SEDONA
Appendix A.4 TARDIS
TARDIS used the Kurucz CD 23 line list to compute the boundbound 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}), metastable 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.
Boundbound atomic data information from the Kurucz CD 1 loneliest for the DDC10 model computed with SEDONA
Boundbound atomic data information from the Kurucz CD 23 line list for all models computed with TARDIS.
References
 Abbott, D. C., & Lucy, L. B. 1985, ApJ, 288, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Abdikamalov, E., Burrows, A., Ott, C. D., et al. 2012, ApJ, 755, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Becker, S. R., & Butler, K. 1992, A&A, 265, 647 [NASA ADS] [Google Scholar]
 Becker, S. R., & Butler, K. 1995a, A&A, 294, 215 [NASA ADS] [Google Scholar]
 Becker, S. R., & Butler, K. 1995b, A&A, 301, 187 [NASA ADS] [Google Scholar]
 Berrington, K. A., Burke, P. G., Dufton, P. L., & Kingston, A. E. 1985, Atomic Data and Nuclear Data Tables, 33, 195 [NASA ADS] [CrossRef] [Google Scholar]
 Blinnikov, S. I. 1996, Astron. Lett., 22, 79 [NASA ADS] [Google Scholar]
 Blinnikov, S. I., & Bartunov, O. S. 1993, A&A, 273, 106 [NASA ADS] [Google Scholar]
 Blinnikov, S. I., Eastman, R., Bartunov, O. S., Popolitov, V. A., & Woosley, S. E. 1998, ApJ, 496, 454 [Google Scholar]
 Blinnikov, S. I., Röpke, F. K., Sorokina, E. I., et al. 2006, A&A, 453, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blondin, S., Dessart, L., Hillier, D. J., & Khokhlov, A. M. 2013, MNRAS, 429, 2127 [NASA ADS] [CrossRef] [Google Scholar]
 Boyle, A., Sim, S. A., Hachinger, S., & Kerzendorf, W. 2017, A&A, 599, A46 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bulla, M., Sim, S. A., & Kromer, M. 2015, MNRAS, 450, 967 [Google Scholar]
 Castor, J. I., Abbott, D. C., & Klein, R. I. 1975, ApJ, 195, 157 [Google Scholar]
 Densmore, J. D., Thompson, K. G., & Urbatsch, T. J. 2012, J. Computat. Phys., 231, 6924 [NASA ADS] [CrossRef] [Google Scholar]
 Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., & Young, P. R. 1997, Astron. Astrophys. Suppl. Ser., 125, 149 [CrossRef] [EDP Sciences] [Google Scholar]
 Dere, K. P., Zanna, G. D., Young, P. R., Landi, E., & Sutherland, R. S. 2019, ApJS, 241, 22 [CrossRef] [Google Scholar]
 Dessart, L., & Hillier, D. J. 2010, MNRAS, 405, 2141 [NASA ADS] [Google Scholar]
 Dessart, L., Hillier, D. J., Blondin, S., & Khokhlov, A. 2014, MNRAS, 441, 3249 [NASA ADS] [CrossRef] [Google Scholar]
 Dessart, L., Audit, E., & Hillier, D. J. 2015, MNRAS, 449, 4304 [Google Scholar]
 Eastman, R. G., & Pinto, P. A. 1993, ApJ, 412, 731 [NASA ADS] [CrossRef] [Google Scholar]
 Ensman, L. M., & Woosley, S. E. 1988, ApJ, 333, 754 [NASA ADS] [CrossRef] [Google Scholar]
 Ertl, T., Woosley, S. E., Sukhbold, T., & Janka, H. T. 2020, ApJ, 890, 51 [CrossRef] [Google Scholar]
 Fleck, J. A. J., & Cummings, J. D. 1971, J. Comput. Phys., 8, 313 [NASA ADS] [CrossRef] [Google Scholar]
 Fontes, C. J., Fryer, C. L., Hungerford, A. L., Wollaeger, R. T., & Korobkin, O. 2020, MNRAS, 493, 4143 [CrossRef] [Google Scholar]
 Hillier, D. J. 1987, ApJS, 63, 947 [NASA ADS] [CrossRef] [Google Scholar]
 Hillier, D. J. 2003, in ASP Conf. Ser., 288, 199 [NASA ADS] [Google Scholar]
 Hillier, D. J. 2012, in From Interacting Binaries to Exoplanets: Essential Modeling Tools, eds. M. T. Richards, & I. Hubeny, 282, 229 [NASA ADS] [Google Scholar]
 Hillier, D. J., & Dessart, L. 2012, MNRAS, 424, 252 [CrossRef] [Google Scholar]
 Hillier, D. J., & Miller, D. L. 1998, ApJ, 496, 407 [NASA ADS] [CrossRef] [Google Scholar]
 Hummer, D. G., Berrington, K. A., Eissner, W., et al. 1993, A&A, 279, 298 [NASA ADS] [Google Scholar]
 Jeffery, D. J. 1999, ArXiv eprints [arXiv:astroph/9907015] [Google Scholar]
 Jerkstrand, A. 2011, PhD thesis, Stockholm University, Sweden [Google Scholar]
 Jerkstrand, A., Fransson, C., & Kozma, C. 2011, A&A, 530, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jerkstrand, A., Fransson, C., Maguire, K., et al. 2012, A&A, 546, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kasen, D., Thomas, R. C., & Nugent, P. 2006, ApJ, 651, 366 [NASA ADS] [CrossRef] [Google Scholar]
 Katz, B., Kushnir, D., & Dong, S. 2013, ArXiv eprints [arXiv:1301.6766] [Google Scholar]
 Kerzendorf, W. E., & Sim, S. A. 2014, MNRAS, 440, 387 [NASA ADS] [CrossRef] [Google Scholar]
 Kerzendorf, W., Sim, S., Vogl, C., et al. 2022, Tardissn/tardis: TARDIS v2022.3.2 [Google Scholar]
 Kingdon, J. B., & Ferland, G. J. 1996, ApJS, 106, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Kozma, C., & Fransson, C. 1992, ApJ, 390, 602 [NASA ADS] [CrossRef] [Google Scholar]
 Kozyreva, A., Shingles, L., Mironov, A., Baklanov, P., & Blinnikov, S. 2020, MNRAS, 499, 4312 [NASA ADS] [CrossRef] [Google Scholar]
 Kramida, A., & Ralchenko, Y. 1999, NIST Atomic Spectra Database, NIST Standard Reference Database 78, type: dataset [Google Scholar]
 Kromer, M., & Sim, S. A. 2009, MNRAS, 398, 1809 [Google Scholar]
 Kurucz, R. L. 2009, in American Institute of Physics Conference Series, 1171, 43 [NASA ADS] [Google Scholar]
 Kurucz, R. L., & Bell, B. 1995, Atomic line list Kurucz CDROM, Cambridge, MA: Smithsonian Astrophysical Observatory, 1995 [Google Scholar]
 Lennon, D. J., & Burke, V. M. 1994, A&AS, 103, 273 [NASA ADS] [Google Scholar]
 Lennon, D. J., Dufton, P. L., Hibbert, A., & Kingston, A. E. 1985, ApJ, 294, 200 [NASA ADS] [CrossRef] [Google Scholar]
 Li, C., Hillier, D. J., & Dessart, L. 2012, MNRAS, 426, 1671 [Google Scholar]
 Liljegren, S., Jerkstrand, A., & Grumer, J. 2020, A&A, 642, A135 [EDP Sciences] [Google Scholar]
 Lucy, L. B. 1999, A&A, 344, 282 [NASA ADS] [Google Scholar]
 Lucy, L. B. 2002, A&A, 384, 725 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B. 2003, A&A, 403, 261 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B. 2005, A&A, 429, 19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lucy, L. B., & Abbott, D. C. 1993, ApJ, 405, 738 [Google Scholar]
 Mazzali, P. A., & Lucy, L. B. 1993, A&A, 279, 447 [NASA ADS] [Google Scholar]
 Mendoza, C. 1983, in IAU Symposium, Planetary Nebulae, ed. D. R. Flower, 103, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Nahar, S. N. 1995, A&A, 293, 967 [NASA ADS] [Google Scholar]
 Nahar, S. N. 2010, NORADAtomicData [Google Scholar]
 Noebauer, U. M., & Sim, S. A. 2019, Liv. Rev. Comput. Astrophys., 5, 1 [CrossRef] [Google Scholar]
 Nomoto, K., Thielemann, F. K., & Yokoi, K. 1984, ApJ, 286, 644 [Google Scholar]
 Nussbaumer, H., & Storey, P. J. 1983, A&A, 126, 75 [NASA ADS] [Google Scholar]
 Nussbaumer, H., & Storey, P. J. 1984, A&As, 56, 293 [NASA ADS] [Google Scholar]
 Pássaro, E. A., E., K. W., A., F., & C., V. 2019, https://doi.org/10.5281/zenodo.4062427 [Google Scholar]
 Pognan, Q., Jerkstrand, A., & Grumer, J. 2022, MNRAS, 510, 3806 [NASA ADS] [CrossRef] [Google Scholar]
 Roth, N., Kasen, D., Guillochon, J., & RamirezRuiz, E. 2016, ApJ, 827, 3 [Google Scholar]
 Seaton, M. J. 1987, J. Phys. B At. Mol. Phys., 20, 6363 [NASA ADS] [CrossRef] [Google Scholar]
 Shen, K. J., Blondin, S., Kasen, D., et al. 2021, ApJ, 909, L18 [NASA ADS] [CrossRef] [Google Scholar]
 Shine, R. A., & Linsky, J. L. 1974, Sol. Phys., 39, 49 [NASA ADS] [CrossRef] [Google Scholar]
 Shingles, L. J., Sim, S. A., Kromer, M., et al. 2020, MNRAS, 492, 2029 [NASA ADS] [CrossRef] [Google Scholar]
 Sim, S. A. 2007, MNRAS, 375, 154 [NASA ADS] [CrossRef] [Google Scholar]
 Spencer, L. V., & Fano, U. 1954, Phys. Rev., 93, 1172 [NASA ADS] [CrossRef] [Google Scholar]
 Sutherland, R. S. 1998, MNRAS, 300, 321 [CrossRef] [Google Scholar]
 Swartz, D. A., Sutherland, P. G., & Harkness, R. P. 1995, ApJ, 446, 766 [NASA ADS] [CrossRef] [Google Scholar]
 Tayal, S. S. 1997a, ApJS, 111, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Tayal, S. S. 1997b, ApJ, 481, 550 [NASA ADS] [CrossRef] [Google Scholar]
 Utrobin, V. P. 2004, Astron. Lett., 30, 293 [NASA ADS] [CrossRef] [Google Scholar]
 van Rossum, D. R. 2012, ApJ, 756, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Verner, D. A., & Yakovlev, D. G. 1995, A&As, 109 [Google Scholar]
 Verner, D. A., Ferland, G. J., Korista, K. T., & Yakovlev, D. G. 1996a, ApJ, 465, 487 [Google Scholar]
 Verner, D. A., Verner, E. M., & Ferland, G. J. 1996b, Atomic Data and Nuclear Data Tables, 64, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Vogl, C., Sim, S. A., Noebauer, U. M., Kerzendorf, W. E., & Hillebrandt, W. 2019, A&A, 621, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Weaver, T. A., Zimmerman, G. B., & Woosley, S. E. 1978, ApJ, 225, 1021 [NASA ADS] [CrossRef] [Google Scholar]
 Wollaeger, R. T., & van Rossum, D. R. 2014, ApJS, 214, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Wollaeger, R. T., van Rossum, D. R., Graziani, C., et al. 2013, ApJS, 209, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Rev. Mod. Phys., 74, 1015 [NASA ADS] [CrossRef] [Google Scholar]
 Wygoda, N., Elbaz, Y., & Katz, B. 2019, MNRAS, 484, 3951 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, H. L., & Pradhan, A. K. 1995a, A&A, 293, 953 [NASA ADS] [Google Scholar]
 Zhang, H. L., & Pradhan, A. K. 1995b, J. Phys. B At. Mol. Phys., 28, 3403 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, H. L., & Pradhan, A. K. 1997, A&AS, 126, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Source code available at https://github.com/artismcrt/artis
CMFGEN, with documentation, is available at www.pitt.edu/hillier
Available at https://github.com/tardissn/carsus
The code is publicly available and can be downloaded from https://www.dropbox.com/sh/kyg1z1xwi8298ru/AAAqzUMbr6AkoVfkSVIYChTLa?dl=8.
CDs 1 and 23 from http://kurucz.harvard.edu/cdroms.html
For the STELLA runs, the files DDC10_29.29s_highres.dat and DDC2 5_62.60s_highres.dat represent a snapshot of the ejecta at ~29 s and ~63 s post explosion, respectively, with 399 shells. The structure is the same but includes four additional columns giving the elemental mass fractions of H, He, N, and P.
The prefactor in parenthesis in the first line of Eq. (10), is equivalent to a simplistic approximation in which a fraction of 0.03 of the energy is emitted in positron kinetic energy which is immediately deposited in the ejecta while the rest is emitted in gammarays with a purely absorptive optical depth given by with a gammaray escape time (e.g. Jeffery 1999) of t_{0} = 40 days.
Data are available online at http://kurucz.harvard.edu
All Tables
Physical ingredients and approximations used in each code for the test models in this paper.
Velocity and temperature at the inner boundary given the requested output luminosity for the TARDIS calculation of the toy06 model.
Model atoms used in CMFGEN calculations up until 10.56 d post explosion for the DDC10 and DDC25 models.
Large Co model atoms used in CMFGEN calculations for the DDC10 (between 11.62d and 15.47 d post explosion) and DDC25 (between 11.62dand 17.02 d post explosion) models.
Model atoms used in CMFGEN calculations from 17.02 d post explosion for the DDC10 model and from 18.72 d post explosion for the DDC25 model.
Large Co model atoms used in CMFGEN calculations for the DDC10 and DDC25 models from 33.15 d post explosion onward.
Reduced Fe IV, Co IV, and Ni IV model atoms used in CMFGEN calculations for the DDC10 and DDC25 models from 40.11 d post explosion onward.
Boundbound atomic data information from the Kurucz CD 1 line list for the toy06 and toy01 models computed with SEDONA
Boundbound atomic data information from the Kurucz CD 1 loneliest for the DDC10 model computed with SEDONA
Boundbound atomic data information from the Kurucz CD 23 line list for all models computed with TARDIS.
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 (^{56}Ni and decay products ^{56}Co and ^{56}Fe, as well as IMEs: Ca S, Si), while for the DDC models only a subset of species are shown for illustration. The ^{56}Ni 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 ^{56}Ni_{0} 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 maximumlight 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 timeweighted 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 Multiband (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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.