StaNdaRT: A repository of standardized test models and outputs for supernova radiative transfer

We present the first results of a comprehensive supernova (SN) radiative-transfer (RT) code-comparison initiative (StaNdaRT), where the emission from the same set of standardized test models is simulated by currently-used RT codes. A total of ten codes have been run on a set of four benchmark ejecta models of Type Ia supernovae. We consider two sub-Chandrasekhar-mass ($M_\mathrm{tot} = 1.0$ M$_\odot$) toy models with analytic density and composition profiles and two Chandrasekhar-mass delayed-detonation 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 in detail the test models, radiative-transfer codes and output formats and provide access to the repository. We present example results of several key diagnostic features.


Introduction
Accurate radiative-transfer (RT) calculations remain a key challenge in the study of astronomical transients such as supernovae (SNe). While advances in computational capabilities and theoretical understanding have allowed great progress in the ability to simulate radiation transport, the large number of physical processes involved, in particular opacity from thousands of atomic transitions with a mixed absorptive and scattering character, prohibit comprehensive 3D calculations based on first principles. Several physical approximations of different forms -in particular different treatments of the significant deviations from local thermodynamic equilibrium (LTE)-are employed by different RT codes to calculate the properties of the gas and of the radiation field. Approximate treatment of atomic physics is also required due to the partially calibrated atomic data.
The back-reaction of radiation on the hydrodynamics provides an additional challenge, requiring the solution of hydrodynamic equations coupled to the RT solution. However, in many cases, and in particular for Type Ia SNe (SNe Ia) at phases beyond several days, which we focus on here, the radiation carries a negligible fraction of the energy and the ejecta are freely expanding homologously. The RT problem in these cases is decoupled from the hydrodynamics problem, the latter providing the initial ejecta profiles ('ejecta models' hereafter). The ejecta pro-files 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, Su-perNu, TARDIS, and URILIGHT; see Sect. 3 for descriptions and references) to create a systematic code-comparison framework. As a first important step, all groups agreed on a set of four test model ejecta and standardised output formats. Each group calculated the resulting radiative display with their respective codes for the same ejecta models and shared the results in a new Article number, page 1 of 25 arXiv:2209.11671v3 [astro-ph.SR] 15 Apr 2023 A&A proofs: manuscript no. main Table 1. Summary of ejecta conditions. The yields for representative species corresponds to the start of the simulations in our model set (2 d post explosion for the toy models and ∼1 d post explosion for the DDC models). The 56 Ni mass is given prior to any decay. 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 sub-Chandrasekhar-mass (sub-M Ch ) toy models with profiles that are defined analytically and two more realistic models that result from hydrodynamical simulations of the M Ch 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.

Test models
The code-comparison test suite consists of four SN Ia models. Two are sub-M Ch models with analytic density and composition profiles (Sect. 2.1; 'toy' models), and the remaining two are M Ch models resulting from hydrodynamical simulations (Sect. 2.2; DDC models). The models were set up or selected based on their 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 low-luminosity 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 d post explosion and initial composition profiles are shown in Figs. 1 and 2, respectively.

Toy models
The toy models were generated using the mk_snia_toy_model.py script (see Sect. 4) using the --highni (for the normal SN Ia model) and --lowni (for the low-luminosity SN Ia model) options. Both models have a total mass M ej = 1.0 M , a kinetic energy E kin = 10 51 erg, and are calculated at the time t f = 2 d post explosion. The models have an exponential density profile (e.g. Jeffery 1999), where v e = E kin 6M ej ≈ 2895 km s −1 (2) is the e-folding velocity, and ρ c = M ej 8πv 3 e t 3 f ≈ 6.32 × 10 −10 g cm −3 is the central density at the chosen time.
The ejecta consist of n spherical shells indexed by i = 1, . . . , n. We use a uniform velocity grid with width ∆v = 50 km s −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 dV = 4πr 2 dr = 4πv 3 e t 3 f z 2 dz, 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: c e −z 0,i (z 2 0,i + 2z 0,i + 2) − e −z 1,i (z 2 1,i + 2z 1,i + 2) z 3 1,i − z 3  2. Composition profiles at the start of our simulations in our model set (2 d post explosion for the toy models and ∼1 d 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.
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 iron-group elements (IGEs), here only 56 Ni and its decay products 56 Co and 56 Fe, and an outer region composed of the intermediate-mass elements (IMEs) Ca, S, and Si with constant mass fractions throughout the layer: X(Ca) = 0.1, X(S) = 0.35, and X(Si) = 0.55. These values were chosen based on the delayed-detonation model DDC10 (see Sect. 2.2 below). Our 'normal' toy model consists of 0.6 M of 56 Ni and 0.4 M of IMEs, while our low-luminosity 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 low-luminosity 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 X IME (m) = 1 − X56 Ni (m). 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 d 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 time-weighted internal energy, tE(t), is equal to the time-integrated time-weighted decay energy deposition rate, tQ(t)dt, as noted by Katz et al. (2013): 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 time-weighted 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 ( 12000 km 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 longer-term radiative display is negligible.

Delayed-detonation models
In addition to the two toy models above, we consider two M Ch delayed-detonation models resulting from hydrodynamical simulations. These were chosen from the DDC model grid presented in Blondin et al. (2013) in order to closely match the 56 Ni yields of the toy models: Our low-luminosity 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 d post explosion by solving the energy equation given by the first law of thermodynamics, assuming the plasma is radiation dominated and neglecting diffusion. Apart from the two-step 56 Ni→ 56 Co→ 56 Fe decay chain, we also treat eight additional two-step 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 one-step 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 d are then remapped onto the 1D, non-LTE, radiative-transfer code CMFGEN of  and evolved until ∼ 1 d post explosion (0.976 d for the DDC10 model and 1.3 d for the DDC25 model; see Sect. 3.2 for details), at which point the ejecta serve as initial conditions for the other radiative-transfer codes presented in Sect. 3.

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

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 three-dimensional and follows the time-evolution of the radiation field and state of the gas. It assumes a strictly homologous velocity law.
Injection of energy into the ejecta is calculated by following the deposition of γ-ray packets that are injected in accordance with the radioactive decays of 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 non-grey and takes into account Compton scattering, photoelectric absorption, and pair-creation opacities. In our standard runs, the code does not include the effects of excitation or ionisation by non-thermal particles. However, Shingles et al. (2020) presented updates to the code that include a Spencer-Fano treatment of non-thermal ionisation as required for late-phase modelling. Results obtained with this improved version (artisnebular) are included for late phases for the models in this study.
The electron temperature in each grid zone is estimated by balancing of heating and cooling rates (accounting for γ-ray and positron deposition, bound-bound, bound-free, and free-free processes). In its regular mode of operation (artis), the code uses an approximate non-LTE treatment to estimate the ionisation state in the ejecta (based on Monte Carlo photoionisation estimators; see Kromer & Sim 2009) and an LTE treatment of excitation is adopted. This approach has been used in most of our published studies, and is also used in most of the artis calculations presented here. However, this method has limitations that become increasingly important at later phases (e.g. it neglects non-thermal heating and ionisation and tends to overestimate the plasma temperature at low densities due to incomplete treatment of cooling by forbidden lines). To improve these issues, Shingles et al. (2020) presented updates to the code that include a full non-LTE population solver (together with the Spencer-Fano solver mentioned above). Results obtained with this improved version (artisnebular) are included for late-phase calculations here.
Monte Carlo estimators are used to track the radiation field in each grid cell. In general, we use volume-based estimators (see Lucy 1999or Noebauer & Sim 2019 to extract radiation-fielddependent 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 black-body radiation field model for bound-bound excitation (see Kromer & Sim 2009 for details). However, the improved artisnebular version uses a more detailed frequency binned representation of the radiation field (see Shingles et al. (2020) for details). The code has the capacity to iterate on each time step with the aim of achieving S. Blondin et al.: Standardised test models and outputs for SN radiative transfer Notes: Column headings: (1) Code name. (2) Numerical method used to solve the radiative-transfer equation: FLD=Flux Limited Diffusion, MC=Monte Carlo, RH-1G = one-group (grey) radiation hydrodynamics, RH-MG = multi-group radiation hydrodynamics, RTE-CMF= Radiation Transfer Equation Co-Moving Frame.
(3) The ejecta are assumed to be in homologous expansion (v = rt) in radiative-transfer codes. This is not the case for radiation-hydrodynamics codes (CRAB, KEPLER, STELLA). (4) Treatment of γ-ray energy deposition. (5) Non-thermal heating, excitation, and ionisation rates are calculated through a solution of the Spencer-Fano equation (Spencer & Fano 1954) or read in from tabulated values (Kozma & Fransson 1992). (6) Solution method for the atomic level populations. LTE(T X ) refers to a solution of the Boltzmann excitation formula setting the temperature to that of the electrons (T e ) or the radiation field (T R ). An approximate non-LTE treatment of excitation scales the Boltzmann occupation numbers by the dilution factor W (cf. dilute-LTE treatment in TARDIS; Sect. 3.9). A non-LTE treatment requires the solution of the rate equations, either including time dependence (dn/dt) or assuming steady-state (statistical equilibrium, dn/dt = 0). (7) Treatment of ionisation. Here LTE(T X ) refers to a solution of the Saha-Boltzmann equation, which can be scaled for an approximate non-LTE treatment (cf. nebular approximation in TARDIS; Sect. 3.9). The non-LTE solution results from the solution of the rate equations, either including time dependence (dn/dt) or assuming steady-state (dn/dt = 0). (8) The radiation field can be computed via a solution of the radiative-transfer equation (possibly assuming steady-state, dJ/dt = 0) or by following the propagation of photon packets in Monte Carlo codes. Alternatively, LTE treatments assume a Planckian radiation field (black body B ν ) at a reference temperature T X , possibly scaled by the dilution factor W. (9) Treatment of line opacity. This can be explicitly line by line, taking into account overlap in the co-moving frame (κ ν ), or with use of the Sobolev approximation. Other treatments involve the use of an approximate frequency-dependent 'expansion' opacity, or assuming a constant value (e.g. KEPLER; Sect. 3.4). (10) Global value of the thermalisation parameter , which sets the probability that a photon absorbed in a given transition is re-emitted in a different transition (see e.g. URILIGHT; Sect. 3.10).
consistency between the radiation field estimates and the packet transport in each step. However, in practise we find that this iteration is not required and we therefore generally simply use the radiation field quantities extracted from the previous time step to estimate the radiative rates that are needed for the current step.
In simulating ultraviolet to infrared photon transport, the code accounts for electron scattering, bound-bound, bound-free, and free-free processes. Bound-free and bound-bound processes are treated using the Macro Atom approach of Lucy (2002Lucy ( , 2003 and adopting the Sobolev approximation for line opacity. The code does not use an expansion opacity (or similar) but treats line opacity based on a frequency-ordered list of transitions treated in the Sobolev limit (i.e. no line overlap is taken into account).
In our simulations, atomic data are primarily drawn from the Kurucz atomic line lists (see Kromer & Sim 2009) -in Appendix A.1 we summarise the ions and numbers of levels and lines that we include. The photon transport is carried out on a 3D Cartesian grid that co-expands with the ejecta. The artis simulations included here were carried out on a 100 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 one to two days.

CMFGEN
CMFGEN is a 1D, non-LTE, time-dependent radiative-transfer code that solves the transfer equation in the co-moving frame of spherical outflows. Details about the code, techniques, and approximations can be found in Hillier & Miller (1998), Hillier (2003, Hillier (2012), and (for SN calculations) in   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, non-local energy deposition from radioactive decay is treated using a Monte-Carlo approach for 3 CMFGEN, with documentation, is available at www.pitt.edu/ hillier γ-ray transport . Non-thermal processes associated with the high-energy electrons produced by Compton scattering and photoelectric absorption of these γ rays are accounted for through a solution of the Spencer-Fano equation (Spencer & Fano 1954;Li et al. 2012).
The temperature structure is constrained through the energy equation that has the form: where D Dt is the Lagrangian derivative, e is the internal energy per unit mass, P is the gas pressure, and˙ decay is the energy absorbed per second per unit volume (see  for further details). We verify the accuracy of the solution by examining a global energy constraint (equivalent to conservation of flux in a static atmosphere; see Hillier & Dessart 2012 for details), and an equation describing energy conservation as applied to the heating and cooling of free electrons. These two equations are related to the above Eq. 9 via the transfer equation, and the rate equations, respectively (e.g. Hillier 2003). Because of model assumptions (e.g. the use of super-levels) these two equations are not satisfied exactly, but the errors (typical at the 1% level or smaller) are too small to affect the validity of the models. Processes contributing to electron heating and cooling include bound-free ionisation and recombination, collisional ionisation and recombination, collisional excitation and de-excitation, freefree emission, Auger ionisation, charge exchange reactions (primarily with H i and He i, and hence negligible in SN Ia ejecta), net cooling from non-thermal processes, and heating by radioactive decay.
Atom and ion-level populations are determined through a solution of the time-dependent rate equations, coupled to the above energy equation and the zeroth and first moments of the radiative-transfer equation (see below). We consider the following processes: bound-bound processes, bound-free processes, collisional ionisation and recombination, collisional excitation and de-excitation, Auger ionisation (Hillier 1987;Hillier & Miller 1998), and non-thermal excitation and ionisation (Li et al. 2012). We additionally consider further processes involving H and He (two-photon decay, charge-exchange reactions, and Rayleigh scattering), although these are negligible here (and completely absent from the toy models, which contain no H or He). To ease the solution of the rate equations, atomic levels are grouped into super levels (see Hillier & Miller 1998 for details).
The frequency-dependent mean intensity J ν is obtained via a solution of the time-dependent 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 so-called Eddington factors f ν = K ν /J ν , where K ν is the second moment of the specific intensity (related to the radiation pressure). The Eddington factors are obtained from a formal solution of the time-independent transfer equation.
We consider the following sources of opacity: electron scattering, bound-free (including photoionisation from excited states), bound-bound 4 , free-free, and Auger ionisation. As noted earlier, Rayleigh scattering and two-photon processes (for H and He only) are also part of the global opacity budget but are negligible here.
A description of the sources of atomic data can be found in the Appendix (Sect. A.2). The number of levels (both superlevels and full levels) and corresponding number of boundbound transitions are given in Tables A.2-A.8. We ignore weak transitions with a g f value 5 below some cutoff, typically set to 10 −4 . For the toy models, the following ions were included: Si iiiv, 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: Ni ii-vii (we also include the ground states of Cl iv, K iii, and V i for the sole purpose of tracking changes in the abundances of radioactive isotopes). For all the aforementioned ions, we also consider ionisations to and recombinations from the ground state of the next ionisation stage (e.g. Ni viii in the case of Ni). As time proceeds and the temperature in the spectrum-formation region drops, the highest ionisation stages have a negligible impact on the RT solution. When this occurs, smaller model atoms are used for these ions, or the ions are removed altogether from the atomic model set.
The toy models were remapped onto a coarser spatial grid with 100 depth points. No remapping was performed for the DDC models. Typical wall-clock runtimes are of the order of 24h per time step on a single computing node with 8-12 CPUs, thus taking 2-3 months to complete for a sequence covering the first 200 d or so post explosion.

CRAB
CRAB is a 1D, implicit, Lagrangian radiation hydrodynamics code developed to study the light curves during SN outbursts and the corresponding outflows evolved from hydrostatic state up to homologous expansion (Utrobin 2004). Non-local energy deposition of γ rays from radioactive decay is determined by solving the corresponding one-group γ-ray transport with the approximation of an effective absorption opacity of 0.06 Y e cm 2 g −1 .
Positrons are assumed to deposit their energy locally. This en-ergy deposition produces non-thermal 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 non-LTE ionisation balance. All elements but H are treated with three ionisation stages. The ionisation balance is controlled by the following elementary processes: photoionisation and radiative recombination, electron ionisation and three-particle recombination, and non-thermal ionisation. In option B, we use an LTE treatment of ionisation and excitation for elements from the option A list or all elements from H to Zn modelled with an arbitrary number of ionisation stages. Atomic and ionic level populations are determined using the Boltzmann formulae and the Saha equations for an element mixture with the local electron or non-equilibrium radiation temperature.
The time-dependent radiation transfer is treated in a 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 limb-darkening effects.
In the inner, optically thick layers, the Rosseland mean opacity is evaluated for the local electron temperature, while in the outer, transparent and semitransparent layers non-LTE effects are taken into account when determining the mean opacities and the thermal emission coefficient. The mean opacities include processes of photoionisation, free-free absorption, Thomson scattering on free electrons, and Rayleigh scattering on neutral hydrogen. The contribution of lines to the expansion opacity is evaluated by the generalised formula of Castor et al. (1975) or by the formula of Blinnikov (1996) using the Sobolev approximation for line opacity. The expansion opacity in an expanding medium is treated with a thermalisation parameter set to 0.9 as recommended by Kozyreva et al. (2020) to model SNe Ia.
The sources of atomic data for processes in continuum can be found in (Utrobin 2004). Oscillator strengths of lines are taken from the Kurucz line database 6 containing nearly 530 000 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 d.

KEPLER
KEPLER is a one-dimensional, implicit Lagrangian hydrodynamics code with appropriate physics for the study of massive stellar evolution and SNe (Weaver et al. 1978;Woosley et al. 2002). The radiation transport is flux-limited diffusion using a single temperature for the matter and radiation.
An important difference between KEPLER and some of the other codes used here is that KEPLER does not assume a coasting configuration. At early times in particular, when the matter is very optically thick, energy input from 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 angle-averaged column depth would be greater than along this radial line. The averaging is therefore approximated by taking a larger opacity. In reality, this number would vary with the radial distribution of 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 flux-limited transport using a single temperature for the radiation and matter, the pressure in the outermost zones can be overestimated. This can lead to a small (of order 10%) overestimate of the conversion of radiation to kinetic energy in those zones at late times. To alleviate this problem, the luminosity is taken to be the maximum of the value at the (electron-scattering) photosphere and the zone that includes 95% of the mass. The former dominated the light curve until after peak.
An important parameter of the calculation is the atomic opacity used to transport thermal radiation. In KEPLER, the total 'optical' opacity consists of two parts: (a) electron scattering, which is calculated using a full Saha solve for all ionisation stages of all 19 species present; and (b) a constant additive opacity, κ a , taken to reflect the effect of Doppler-broadened lines. The 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 zero-age 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 odd-Z 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. High-velocity zones would also cause difficulty with editing the luminosity at late times if the light-crossing time ceased to be negligible. Removing this high-velocity material decreased the kinetic energy of all models by about 2%. This should have negligible effects on the light curve. The zoning of the DDC models was relatively coarse by traditional KEPLER standards. After trimming the high-velocity zones, only 80 zones remained. The zoning of the toy models was better with roughly 600 zones remaining. Rezoning was not carried out. For a light curve with no nuclear burning, the zoning is adequate. All calculations took at most a few minutes on a laptop.

SEDONA
SEDONA is a time-dependent, multi-frequency Monte Carlo radiative transfer code originally developed to study SN light curves, spectra, and polarisation (Kasen et al. 2006). For this comparison study, the gas properties were tracked on 1D spherical Lagrangian zones and are evolved under the assumption of homologous expansion.
The radioactive decays of 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 radioactive-decay processes were used in this comparison study.
The γ-ray packets were transported using the Monte Carlo procedure subject to simplified treatments of bound-free absorption and Compton down-scattering. The γ-ray interactions were treated as heating terms that entered into the thermal heating balance that sourced a population of thermal photon packets.
The temperature in each shell was computed assuming radiative equilibrium. In more detail, the LTE assumption let us set the gas emissivity to the Planck function multiplied by the frequency-dependent line opacity. The temperature was then iteratively adjusted until the frequency-integrated emissivity was equal to the total radiative energy from both γ-rays and reemitted photon packets absorbed over the previous particle propagation step, as measured in the co-moving frame of the fluid.
For most of its published applications, SEDONA assumes LTE in order to compute opacities, although non-LTE capabilities have been developed and were recently applied (Roth et al. 2016;Shen et al. 2021). For this comparison study, only the LTE opacity mode was used in order to compute the excitation and ionisation states of the gas. More precisely, this means that ionisation fractions for each element were computed by simultaneously solving the Saha equation using a local temperature, and the charge conservation equation across all elements and isotopes present. Meanwhile, the bound-electron level populations within each ionisation stage were set by the Boltzmann factors given by the local temperature.
Photon packets were transported in three dimensions as a direct solution to the time-dependent radiative transfer equation. For all interactions with the gas, the packets were mapped to the 1D Lagrangian zones.
Once the level populations had been computed, the 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 bound-bound opacities and electron scattering opacity, using the Thomson cross section. For DDC10, simplified bound-free and free-free opacities were also included in addition to the previously mentioned opacities, although they did not have a noticeable effect during the early stages of the explosion, which is what we wish to compare using these calculations.
A thermalisation parameter of 1.0 was used for the toy01 and toy06 models. This means that all photon packets (other than γ-ray packets) that were absorbed were immediately thermalised, so that each absorption event was followed by 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 re-sampled.
The atomic data used for the bound-bound transitions were taken from Kurucz CD 1. This is a larger set of atomic data than CD 23. The details are described in Appendix A.3.
For this comparison study, the co-moving frequency grid for the thermalised photons (i.e. not the γ-rays) used 17664 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 post-explosion required about 10 CPU hours, and could be efficiently performed in parallel across several hundred processors, reducing the elapsed wall time to less than one hour.

STELLA
The multi-group radiation hydrodynamics code STELLA (Blinnikov & Bartunov 1993;Blinnikov et al. 1998Blinnikov et al. , 2006 is capable of computing the evolution of the radiation field coupled to the hydrodynamics, as well as the bolometric light curve, spectral energy distribution, and resulting broad-band magnitudes and colours. Therefore, STELLA does not require the condition of homologous expansion and is capable of treating shock propagation and any dynamical processes taking place in the SN ejecta. Energy deposition from 56 Ni and 56 Co radioactive decay is treated in a one-group 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 least-squares method.
The opacity includes photoionisation, free-free 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 lower-resolution grid (with 202 zones), while the DDC models were computed at a higher resolution (with 399 zones). The typical runtime on a single processor is a few hours at most.

SUMO
SUMO (Jerkstrand 2011;Jerkstrand et al. 2011Jerkstrand et al. , 2012) is a homologous non-LTE code with radiative transfer. It is specialised to calculate spectra and light curves in the post-peak phases of the SN. The code is written in Fortran and is parallelised with MPI.
Gamma-ray transfer is done by ray tracing using a grey opacity of 0.06Y e cm 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 non-thermal electrons following the scattering of gamma rays and positrons is solved for with the Spencer-Fano method (Kozma & Fransson 1992).
Zone temperatures are solved from the first law of thermodynamics considering heating by non-thermal processes and photoionisation, and cooling by net collisional line excitation, freefree emission, and recombination. The temperature is solved either in steady state (heating = cooling) (Jerkstrand et al. 2012) or time-dependently (Pognan et al. 2022). Here, the steady-state mode was used.
The rate equations are solved by considering spontaneous radiative decay with Sobolev escape probabilities (assuming homology) modified to include continuum processes and strong line overlaps, thermal and non-thermal collisions (excitations and ionisations), photoexcitation and de-excitation, photoionisation, and recombination. Options exist to also include charge transfer processes and molecular chemistry (Liljegren et al. 2020), but these were not considered here. The rate equations can be solved either in steady state (inflow = outflow) (Jerkstrand et al. 2012) or time dependently (Pognan et al. 2022). Here, the steady-state mode was used.
Radiative transfer is calculated with a Monte Carlo simulation that is iterated with the solvers for temperature, ionisation, and excitation. During the transfer, photon packets can experience electron scattering, bound-free and free-free absorption, and line absorption. The line transfer is resolved line by line (with Sobolev formalism) rather than using an expansion opacity formalism. Photoexcitation is either fully coupled to the non-LTE solutions (for low-and mid-lying levels), or decoupled (for high-lying levels), instead giving a fluorescence cascade on the spot. The radiation field is computed in steady state.
The atomic data come from a variety of sources, mostly described in Jerkstrand et al. (2011) andJerkstrand et al. (2012). Model atoms have LS-states resolved, that is they do not use super levels. In models computed up to 2022 (including the ones here), ions up to and including doubly ionised were included, although higher ions are now being included (Pognan et al. 2022). For the toy models computed here, the model atoms are Fe i (496 levels), Fe ii (578 levels For the runs here, the ejecta were resampled to ∆v = 500 km s −1 , and truncated at v = 30000 km s −1 . The radiative transfer was computed on a wavelength grid from 400 Å to 25000 Å, with a logarithmic resolution dλ/λ = 10 −3 . Typical wall-clock runtimes for a single epoch are a few hours on a typical 128 core setup.

SuperNu
SuperNu is a multi-group LTE radiative-transfer code that employs Implicit Monte Carlo (IMC) and Discrete Diffusion Monte Carlo (DDMC) (Wollaeger et al. 2013;Wollaeger & van Rossum 2014). IMC solves the thermal radiative-transfer equations 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 mean-freepath 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 semi-analytic radiative transfer tests (Wollaeger et al. 2013) and on the W7 model of SNe Ia (Nomoto et al. 1984;Wollaeger & van Rossum 2014).
For the γ-ray transfer, SuperNu employs a constant absorption opacity of 0.06 Y e cm 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 semi-implicit linearisation (Fleck & Cummings 1971) of the comoving internal energy equation (Wollaeger et al. 2013): internal energy is recast to gas temperature using the standard relation ∂e = c v ∂T , where e is internal energy, c v is heat capacity at constant volume, and T is gas temperature. Due to the IMC time linearisation, energy deposited from gamma-rays and (locally) from beta particles appears simultaneously in both the comoving internal energy equation and the radiation equation.
Ionisation and excitation are both treated with Saha-Boltzmann statistics evaluated at the gas density and temperature at the beginning of the time step. The multi-element system is solved iteratively by converging the free electron number. The resulting population densities are then used to calculate opacities. The radiation field is represented by fully time-dependent Monte Carlo packets, which are sourced from the LTE emissivity. The radiation field is not constrained to be in equilibrium with the gas, and so in general the system is 'two-temperature'.
Opacity is discretised into groups via direct integration over co-moving wavelength (but see Fontes et al. 2020 for a study with SuperNu using expansion opacity). Opacity in SuperNu includes free-free (Sutherland 1998), and bound-free (Verner & Yakovlev 1995;Verner et al. 1996a) processes, as well as the bound-bound opacities from the Kurucz line lists 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 32000 Å.
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.

TARDIS
TARDIS is an open-source steady-state 1D radiative-transfer code that uses indivisible energy packets as its transport quanta following the methods in Abbott & Lucy (1985), Lucy & Abbott (1993), Mazzali & Lucy (1993), Lucy (2002), Lucy (2003), and Lucy (2005). Kerzendorf & Sim (2014) describes the initial version of the code which was primarily used to model SNe Ia. Subsequently, the code has been significantly enhanced to include a non-thermal approximation treatment for helium (Boyle et al. 2017) and the continuum processes and relativity treatments required for hydrogen-rich SNe (Vogl et al. 2019). TARDIS has been continuously enhanced since then (see e.g. Kerzendorf et al. 2022, similarly gamma-ray energy deposition is a new module of TARDIS, but has not been used for the models in this work). Full documentation and an extended physics walkthrough for TARDIS can be found online 9 .
TARDIS assumes a steady-state homologously expanding SN envelope and injects Monte Carlo packets -randomly sampled from a given distribution (by default a black-body)-from an inner boundary. The code supports bound-bound, boundfree, free-free, and Thomson opacities with several redistribution schemes from simple scattering to a macro-atom (Lucy 2002(Lucy , 2003. Summary packet statistics are used to estimate radiation field quantities (temperature, dilution factors, mean intensities, heating-and photo-ionisation rates, and line source functions). The estimated quantities are used to calculate the ionisation and excitation populations in steady state with a choice between A&A proofs: manuscript no. main Table 3. Velocity and temperature at the inner boundary given the requested output luminosity for the TARDIS calculation of the toy06 model. Notes: L req = requested output luminosity; IB = inner boundary.
LTE and several formulations of non-LTE (nebular approximation of Abbott & Lucy 1985 for ionisation, and dilute-LTE and full non-LTE for excitation). TARDIS then calculates Sobolev optical depths for line interaction, and opacities for continuum processes. These values are then used in the subsequent Monte Carlo step, which produces summary statistics for updating the opacities. The other convergence criterion is a match of the integrated packet output luminosity and the requested luminosity, which is achieved through iterative adjustments of the inner temperature.
The setup files for TARDIS that are used in this work are available online 11 using an atomic data set from Kurucz CD 23. In the setup used for the code comparison work, we run TARDIS in a mode that self-consistently finds a temperature stratification given an inner boundary velocity and output luminosity. We used the mean bolometric luminosity of the other comparison codes for our output luminosity. The structure of the toy01 and toy06 models was set at 520 shells between 9000 and 35000 km s −1 . The DDC models were truncated to 40 shells between 9000 and 35000 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 dilute-lte 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 20000 Å with 2 Å bin width. The models were run on one CPU with runtimes of less than one hour.

URILIGHT
URILIGHT is a time-dependent Monte-Carlo code written in Fortran 90 by Yoni Elbaz based on the approximations that are used in SEDONA (see Sect. 3.5 above, Kasen et al. 2006 and references therein), in particular assuming homologous expansion. A detailed description of this program and previous comparisons to other published radiativetransfer codes for several benchmark problems are presented in Wygoda et al. (2019). The code is publicly available and can be downloaded from https://www.dropbox.com/sh/ kyg1z1xwi0298ru/AAAqzUMbr6AkoVfkSVIYChTLa?dl=0.
Energy deposition resulting from the decay of radioactive isotopes is calculated by a Monte-Carlo solution of the γ-ray transport, for which interaction with matter is included through Compton scattering and photoelectric absorption. For the calculations in this work, only 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 Boltzmann-distribution.
Opacities include bound-bound and free-free absorption and Thomson scattering off free electrons. The atomic line data for the bound-bound transitions, which constitutes the main and most important source of opacity in SNe Ia, are taken from the extended set of lines by Kurucz 12 . Following a bound-bound interaction, most photons (a fraction , which is a global parameter of the simulation; see Kasen et al. 2006) are thermalised and reemitted at a different wavelength with a distribution set by the emissivity. The rest of the photons (fraction 1 − ) are re-emitted at the same wavelength (within the line width). As in Wygoda et al. (2019), we used = 0.8 in the runs performed here.
The runs here use 162 spatial cells for the toy01 and toy06 models and 115 spatial cells for the DDC10 and DDC25 models. All models are run with 128 time steps logarithmically spaced between 2 and 210 days, and a uniform spectrum resolution of 10 Å between 100 and 30000 Å. Each run typically took 30 hours on a single core, or approximately one hour when parallelised.

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 https://github. com/sn-rad-trans/data1. 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).

Ejecta model files
As described in Sect. 2, the code comparison is performed using four SN Ia models (main parameters in Table 1). The RT input files that were distributed among the groups are provided in the repository, including two toy-model files snia_toy01_2d.dat and snia_toy06_2d.dat and the two delayed-detonation model files DDC10_0.976d.dat and DDC25_1.300d.dat. Notes: RT Method gives the numerical method used to solve the radiative-transfer equation (see Table 2).

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

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.

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: Each block content consists of NVEL rows with the following NSTAGES+1 columns (7 in this example): 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

Useful Python codes
All the figures presented in this paper are automatically generated and accessible in a Python notebook all_plots.ipynb which is available in the data repository. A separate Python notebook photometry.ipynb in the code-comparison1/ directory is used to generate the synthetic photometry (wsynphot_mags files) from the spectra files on the fly. Three further useful Python codes can be used to: 1. Generate the toy models Location: data1/input_models/mk_snia_toy_model.py Calling syntax: python mk_snia_toy_model.py --highni python mk_snia_toy_model.py --lowni This code generates the toy models used in this paper, as explained in Sect. 2.1. The upper command generates the toy06 model (snia_toy06_2d.dat), while the lower line generates the toy01 model (snia_toy01_2d.dat). Comments inside the code provide information on how to create new similar test models with different parameters.

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.

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).

Example results
In this section we provide example results that are extracted from the outputs of the different simulations. The purpose of this section is to illustrate the contents of the data repository and comparisons that can be made using it. No comparison to observations or in-depth investigations of the sources of differences in results are made here. We note that while some groups provided a few sets of results for different physical approximations that appear in the repository, only one set of results is shown for each code in this set of examples. One exception is ARTIS, for which both the regular (artis) and nebular (artisnebular) versions are included for each model in order to present results for this code from early to late times. The results are provided in the following subsections. In Sect. 5.1, bolometric properties of the emission are shown for simulations of the toy06 model, including the pseudobolometric (UVOIR) light curve and the energy deposition rate as a function of time. Rise times and peak luminosities are provided for the toy06 and toy01 models. In Sect. 5.2, the resulting profiles of temperature and mean ionisation state of cobalt are shown at selected times. In Sect. 5.3, multi-band light curves and colour curves are shown for model toy06. B-band rise times are provided for all models where available. Example spectra at selected  5. Total energy deposition rate from γ rays and positrons for model toy06, normalised to the analytic function given by Eq. 10 (the normalisation allows differences to be seen more clearly). 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.

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   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: .45 e −t/8.8d + 1.45 e −t/111.3d 10 43 erg s −1 where M Ni = 0.6 M is the total (undecayed) 56 Ni mass in the toy06 model 15 . 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): 15 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 gamma-rays with a purely absorptive optical depth given by τ γ = t 2 0 /t 2 with a gamma-ray escape time (e.g. Jeffery 1999) of t 0 = 40d. the instantaneous ratio of the energy deposition rate and the luminosity (left) and the ratio of cumulative time-weighted integrals of the energy deposition rate and luminosity (right).

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=0 i · f i , 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).

Multi-band light curves and colours
In this section, multi-band properties are shown for the different codes for different models. The photometry is extracted from the spectra reported in the spectra output files using the wsynphot package 16 (using the Vega calibration spectrum alpha_lyr_stis_003.fits 17 and a third-order spline interpolation). In Fig. 8, the U BVRI JHK light curves are shown for model toy06 up to 100 days after explosion. The B-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.

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-10000 Å to allow direct comparison. This does not represent the full range available in the output files, which differs between different codes. In Fig. 11, early-time spectra (5 days post explosion and around peak luminosity) are shown while in Fig. 12 we show spectra at later times (50 days and 200 days post explosion). We note that nebular times require specialised treatment in non-LTE which is only implemented in the ARTIS, CMFGEN, and SUMO codes for this paper 18 .

Discussion and outlook
In this paper, we present an online public data repository that includes test models for radiative-transfer (RT) calculations of emission from SNe and standardised simulation outputs by ten different groups that allows direct comparison. Python scripts that generate the analytic toy models and read the different formats are provided as well.
The main purpose of the repository is to allow studies of the different physical approximations involved in the different codes and to assess the robustness of different predictions of radiation transfer. In addition, the repository can be used for finding bugs in the codes or can provide checks in the development of new codes.
We plan to extend the set of test models to include other SN types, as well as multi-dimensional models. In addition, we in-16 https://github.com/starkit/wsynphot 17 https://archive.stsci.edu/hlsps/reference-atlases/ cdbs/calspec/alpha_lyr_stis_003.fits 18 SEDONA also has non-LTE capabilities but these were not used here. tend 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.    Fig. 11 and is not repeated here for sake of clarity.
Tables A. 2-A.8 give the number of levels (both super-levels and full levels; see Miller 1998 andHillier 2010 for details) for the model atoms used in the radiative-transfer calculations presented in this paper. N SL refers to the number of super levels used for the solution of the rate equations, and N full refer the number of full levels used to solve the transfer equation and compute the observed spectrum. We report the uppermost level for each ion treated in the fourth column. 'W' refers to states in which higher states (usually f or higher) have been combined into a single level. In the last column we give the number of bound-bound transitions in the model ion taking into account all N full levels. Ions for which N SL = N full = 1 (with no entries for the last level configuration or number of lines) correspond to the final ionisation stage of a given element, for which ionisations to and recombinations from the ground state are considered. The ions Cl iv, K iii, and V i were included for the sole purpose of tracking changes in abundance of radioactive isotopes. The entries for those ions in Tables A. 2-A.4 also have N SL = N full = 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. These authors also provide transition probabilities to states in the ion continuum. The largest source of oscillator data is from Kurucz (2009) 19 ; its principal advantage over many other sources (e.g. Opacity Project) is that LS coupling is not assumed. More recently, non-LS oscillator strengths have become available through the Iron Project (Hummer et al. 1993), and work done by the atomic-data group at Ohio State University (Nahar 2010). Other important sources of radiative data for Fe include Becker & Butler (1992; Nahar (1995). Energy levels have generally been obtained from the National Institute of Standards and Technology. Collisional data are sparse, particularly for states far from the ground state. The principal source for collisional data among low-lying states for a variety of species is the tabulation by Mendoza (1983); other sources include Berrington et al. (1985), Lennon et al. (1985), Lennon & Burke (1994), Shine & Linsky (1974), Tayal (1997a,b), Zhang & Pradhan (1995a,b, 1997. Photoionisation data is taken from the Opacity Project (Seaton 1987) and the Iron Project (Hummer et al. 1993). Unfortunately Ni and Co photoionisation data are generally unavailable, and we used crude approximations. Charge exchange cross-sections are from the tabulation by Kingdon & Ferland (1996).  Notes: Due to a g f cut (level dependent, g f > 10 −4 ) only 377 383 lines were included in the non-LTE calculations of the level populations. 692 911 lines were include when computing the observed spectrum. Prior to 2.4 d post explosion a higher cut was used (g f > 10 −3 ) in order to ease convergence. From 25.93 d onward we omitted the highest ionisation stages Fe v, Co v, and Ni v since their populations are too low to affect the radiative transfer.  Notes: We only report differences with respect to the model atoms shown in Table A.4 (the total values in the last line take into account all ions). In these calculations, we omitted the ions C iv, O iv, Fe vi-vii, Co vi-vii, and Ni vi-vii since their populations were too low to affect the radiative transfer. Due to a g f cut (level dependent, g f > 2 × 10 −3 ) only 206 738 lines were included in the non-LTE calculations of the level populations. 432 953 lines were include when computing the observed spectrum.  Notes: We only report differences with respect to the model atoms shown in Table A.6 (the total values in the last line take into account all ions). Due to a g f cut (level dependent, g f > 10 −4 ), only 1 277 321 lines were included in the non-LTE calculations of the level populations. We note that 1 860 140 lines were included when computing the observed spectrum. Notes: We only report differences with respect to the model atoms shown in Table A.6 (the total values in the last line take into account all ions). Due to a g f cut (level dependent, g f > 10 −4 ), only 1 118 403 lines were included in the non-LTE calculations of the level populations. We note that 1 663 922 lines were include when computing the observed spectrum.  Table A.11 tabulates the total number of levels (N levels ), meta-stable levels (N meta ), and lines (N lines ) for the atoms used in the four models. Atoms present in the atomic data file but not used by the models were not listed.
Article number, page 23 of 25 A&A proofs: manuscript no. main