Spectral modeling of type II supernovae
I. Dilution factors
^{1}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741
Garching, Germany
email: cvogl@mpagarching.mpg.de
^{2}
Astrophysics Research Centre, School of Mathematics and Physics, Queen’s University Belfast, Belfast, BT7 1NN
UK
email: s.sim@qub.ac.uk
^{3}
European Southern Observatory, KarlSchwarzschildStr. 2, 85741
Garching, Germany
^{4}
Physik Department, Technische Universität München, JamesFranckStr. 1, 85741
Garching, Germany
Received:
22
June
2018
Accepted:
1
November
2018
We present substantial extensions to the Monte Carlo radiative transfer code TARDIS to perform spectral synthesis for type II supernovae. By incorporating a nonLTE ionization and excitation treatment for hydrogen, a full account of free–free and boundfree processes, a selfconsistent determination of the thermal state and by improving the handling of relativistic effects, the improved code version includes the necessary physics to perform spectral synthesis for type II supernovae to high precision as required for the reliable inference of supernova properties. We demonstrate the capabilities of the extended version of TARDIS by calculating synthetic spectra for the prototypical type II supernova SN1999em and by deriving a new and independent set of dilution factors for the expanding photosphere method. We have investigated in detail the dependence of the dilution factors on photospheric properties and, for the first time, on changes in metallicity. We also compare our results with the previously published sets of dilution factors and discuss the potential sources of the discrepancies between studies.
Key words: radiative transfer / methods: numerical / stars: distances / supernovae: general / supernovae: individual: SN1999em
© ESO 2019
1. Introduction
In recent years the availability of spectral data for hydrogenrich supernovae (Type II; SNe II) has increased dramatically. Measurements for hundreds of SNe II are now publicly accessible (see, e.g., Poznanski et al. 2009; D’Andrea et al. 2010; Hicken et al. 2017, for recent data releases), providing a dataset that contains a wealth of information about the kinematics of the explosion, the progenitor systems (e.g., Jerkstrand et al. 2012), the circumstellar material (e.g., Quimby et al. 2007), and much more. Most of the analysis of this data has focused on the study of easily measurable spectral parameters such as line absorption velocities or equivalent widths and their correlations (see, e.g., Gutiérrez et al. 2017a,b), omitting the full information contained in the spectra. To establish connections between these parameters and the underlying quantities, such as the metallicity, most studies rely on approximate relations that have been calibrated based on theoretical models (see, e.g., Anderson et al. 2016). In contrast, only a few wellobserved type II supernovae, such as SN1999em (Baron et al. 2004; Dessart & Hillier 2006), SN2005cs (Baron et al. 2007; Dessart et al. 2008) or SN2006bp (Dessart et al. 2008), have been studied using detailed radiative transfer models, which provide a direct way to infer information about the chemical composition, the density profile and other parameters from the full spectral time series. This applies in particular to the use of SNe II as distance indicators, despite the fact that an absolute distance estimate is a natural byproduct of a quantitative spectroscopic analysis (see Baron et al. 1995, 1996a, 2004, 2007; Lentz et al. 2001; Mitchell et al. 2002). SNe II have a long history as cosmological probes (Kirshner & Kwan 1974; Schmidt et al. 1994) and have regained popularity in recent years due to the increased availability of data at high redshifts (e.g., Poznanski et al. 2009; Gall et al. 2018) and, in the era of high precision cosmology, due to the rising need for independent tests of our cosmological models. Recent efforts have focused mainly on methods that rely on various observed correlations between photometric and spectroscopic parameters such as the standard candle method (SCM) by Hamuy & Pinto (2002), the photometric color method by de Jaeger et al. (2015) or the photospheric magnitude method by Rodriguez et al. (2014). Both Gall et al. (2018) and de Jaeger et al. (2017) demonstrate that with the SCM or the expanding photosphere method (EPM) distance measurements of SNe II up to redshifts of ≈0.34 are feasible, highlighting the progress that has been made possible through the availability of new data. In contrast, the determination of distances from radiative transfer modeling has stagnated in recent years. Both the tailored EPM (Dessart & Hillier 2006; Dessart et al. 2008) and the spectral fitting expanding atmosphere method (SEAM) (Baron et al. 1995, 1996a, 2004, 2007) have never been applied outside the local universe. Nevertheless, their independence of the cosmic distance ladder as well as their foundation in wellunderstood physics make them a promising independent tool for cosmology.
Motivated by the wealth of available spectral data and the unique diagnostic abilities of radiative transfer modeling, we have developed a new numerical tool for performing spectroscopic analysis of SNe II. Since our main goal is to provide a tool for parameter inference, we neglect timedependent effects in favor of computational expediency. Currently, the high computational costs prevent the application of numerical methods that selfconsistently simulate the time evolution of the radiation field and the plasma state based on initial conditions (Dessart & Hillier 2011; Dessart et al. 2013) to this purpose. Our approach is an extension of the Monte Carlo radiative transfer code TARDIS (Kerzendorf & Sim 2014), which was originally developed for spectral synthesis in type Ia supernovae (SNe Ia). We have extensively modified and improved the physical treatment of radiative transfer implemented in TARDIS to be applicable to the modeling of SNe II atmospheres. This improved version of the code is then used to calculate a new, independent set of dilution factors for the EPM. In the EPM the dilution factors as introduced by Hershkowitz et al. (1986a,b) and Hershkowitz & Wagoner (1987) correct for the deviation of the supernova emission from that of a blackbody of the same color temperature. They provide the possibility to compare our model calculations to previously published numerical results by (Eastman et al. 1996; E96 from now on and Dessart & Hillier 2005a; D05 from now on) in a simple parametrized fashion. Currently, the systematic discrepancies between the two sets of dilution factors constitute one of the most significant sources of uncertainty in the EPM, accounting for differences of roughly 20% in the inferred distance (e.g., Takats & Vinko 2006; Jones et al. 2009; Gall et al. 2016, 2018). This significant uncertainty highlights the need for additional calculations based on independent numerical methods to understand and resolve the current tension.
The structure of the paper is as follows. We begin with a detailed description of the physical extensions and their numerical implementation into TARDIS in Sect. 2. In Sect. 3, we provide a brief review of the EPM and discuss the basic physics of the dilution factors. As a first application of the extended version of TARDIS, we present radiative transfer models for two epochs of the prototypical SN II SN1999em in Sect. 4. The next sections are dedicated to the presentation and discussion of our main application, the calculation of a new set of EPM dilution factors. The setup of the necessary grid of supernova models is described in Sect. 5, followed by an analysis of the calculated dilution factors in Sect. 6. Here, we focus particularly on the differential influence of the model parameters such as photospheric density or metallicity. To put our results into context and to understand the differences between the published set of dilution factors, a comparison to previous studies is given in Sect. 7. We investigate the differences in the adopted numerical approaches and examine the different choices for the atmospheric properties. Finally, we summarize our results and give an outlook in Sect. 8.
2. Method
We present an extended version of the onedimensional Monte Carlo (MC) radiative transfer code TARDIS (Kerzendorf & Sim 2014) that has been significantly extended for the application to SNe II. TARDIS is based on the indivisible energy packet MC methods of Lucy (1999a,b, 2002, 2003) and has been developed for rapid spectral modeling of SNe Ia. It has been used to study various aspects of SN Ia explosion physics. Applications include abundance tomographies (Barna et al. 2017), a study of spectral signatures of helium in doubledetonation models (Boyle et al. 2017), as well as analyses of SNe Iax spectra (Magee et al. 2016, 2017). In these studies only the effects of Thomson scattering and bound–bound line interactions are simulated in detail. This is a reasonable approximation for SNe Ia but not for SNe II, which have a higher ratio of continuum to line opacity due to the hydrogenrich composition. To adapt TARDIS to these conditions, we extend our treatment of radiation–matter interactions to include boundfree, free–free as well as collisional processes using the macro atom scheme of Lucy (2002, 2003) as outlined in Sect. 2.1. Further necessary improvements to the code can be motivated based on the peculiarities of radiative transfer in SNe II. SNe II atmospheres are characterized by comparatively low densities at the photosphere and a scattering dominated opacity. Due to the low densities, collisions are ineffective at coupling the level populations and ionization and excitation are mainly controlled by the radiation field. The radiation field is dilute compared to its equilibrium value as a result of the dominance of electronscattering opacity and thus significant departures from local thermodynamic equilibrium (LTE) arise even far below the photosphere (see, e.g., Dessart & Hillier 2011). To address this issue, we have extended the code as outlined in Sect. 2.2. Another consequence of the scattering dominated environment is that relatively high optical depths on the order of τ ∝ 𝒪(10) are needed to guarantee a full thermalization of the radiation (see, e.g., Eastman et al. 1996). At such high optical depths the atmospheric structure is strongly affected by relativistic transfer effects as demonstrated by Hauschildt et al. (1991). The inclusion of these effect in TARDIS is described in Sect. 2.1.4.
2.1. Monte Carlo simulations
To find a consistent solution for the plasma state and the radiation field, TARDIS performs a series of Monte Carlo radiative transfer simulations. At every radiative transfer step, a large ensemble of indivisible energy packets (see Abbott & Lucy 1985; Lucy 1999a, 2002, 2003) representing monochromatic photon bundles is initialized at the inner boundary. Initial packet properties are assigned under the assumptions of the LTE diffusion limit. Thus, packet frequencies are sampled from a blackbody distribution at the inner boundary temperature T_{i} and propagation directions are selected according to zero limbdarkening in the comoving frame. Uniform packet energies are chosen such that the injected packets carry a total comoving frame luminosity , where σ is the StefanBoltzmann constant and R_{i} is the radius of the inner boundary. With initial properties assigned, the propagation of the packets is simulated under the assumption of a steadystate, that is to say, neglecting time dependence, as outlined in the following section.
2.1.1. Packet propagation
After initialization, each packet is followed until it leaves the computational domain through the inner or outer boundary. Between the boundaries the supernova atmosphere has been discretized into equidistant, spherical shells. Within each shell the plasma properties such as the opacity or the electron temperature are assumed to be constant. During the propagation the effects of Thomson scattering, hydrogen boundfree, free–free, bound–bound as well as collisional processes are taken into account. As described in Kerzendorf & Sim (2014), line opacity is treated in the Sobolev approximation (see Sobolev 1957). For microturbulent velocities on the order of 100 km s^{−1}, this is as accurate as the comoving frame method in describing the formation of the Balmer lines in SNe II (see Duschinger et al. 1995). Following Lucy (2003), the free–free opacity
is evaluated with free–free gaunt factors set to unity. Here, N_{j, k} denotes the number density of ionization stage j of element k, T_{e} is the electron temperature, n_{e} the electron density and ν the frequency. The prefactor α_{ff} has the value 3.69 × 10^{8} cm^{5} K^{1/2} s^{−3}. Since hydrogen is the dominant source of boundfree opacity in SNe II, we restrict the inclusion of these processes currently only to this element. However, since an extension to more species is conceptually straightforward, we present the governing equations in their general form. Thus, the opacity resulting from photoionizations of electrons in level i of ion j, k is given by
where n_{i, j, k} and denote the actual and the respective LTE level number densities (see Eq. (5.25) of Hubeny & Mihalas 2014). The crosssection for photoionzation α_{i, j, k → j + 1, k}(ν) is obtained from tabulated values through linear interpolation.
To account for the inclusion of hydrogen boundfree, free–free as well as collisional processes small modifications to the packet propagation procedure of Kerzendorf & Sim (2014, Sect. 2.6) have been necessary. In particular for continuum interactions an additional MC experiment is needed to determine the physical absorption mechanism. The probabilities for bound free absorption, free–free absorption and Thomson scattering are given by χ^{bf}/(χ^{bf} + χ^{ff} + χ^{Th}), χ^{ff}/(χ^{bf} + χ^{ff} + χ^{Th}) and χ^{Th}/(χ^{bf} + χ^{ff} + χ^{Th}) respectively. If a boundfree process is selected, a specific continuum for absorption has to be assigned according to the probabilities for photoionization from specific levels i of ion j, k. Regardless of the type of interaction, we use the macro atom scheme of Lucy (2002, 2003) to select an emission channel as outlined in Sect. 2.1.2. For boundfree and free–free emission, the packet has to be assigned an appropriate frequency before the propagation can be resumed. We employ the approximate sampling rule of Lucy (2003, Eq. (41)) for free–free processes and linear interpolation on precomputed values of the emissivity for boundfree interactions.
2.1.2. Macro atom
We use the macro atom scheme of Lucy (2002, 2003) for a general treatment of complicated radiationmatter interactions, such as recombination cascades, fluorescent line emission or cooling emission. In this scheme, packet splitting for processes with multiple emission channels is avoided by assigning the total energy of the packet randomly to one possible interaction channel according to a set of rules derived from the assumption of statistical equilibrium. In Kerzendorf & Sim (2014), only the redistribution of excitation energy created by bound–bound absorption events was simulated using the macro atom machinery. We introduce indivisible packets of thermal kinetic energy (kpackets) and ionization energy (ipackets) in addition to the packets of excitation energy (macro atoms) included in Kerzendorf & Sim (2014) to treat continuum interactions. kpackets can be created by boundfree and free–free absorption events as well as collisional deactivations of ipackets or macro atoms. Since both thermal and ionization energy are created in photoabsorption events, the rpacket is transformed into a kpacket with probability p^{k} = ν_{i, j, k}/ν′ and into an ipacket otherwise. Here, ν_{i, j, k} is the threshold for ionization and ν′ is the frequency of the rpacket in the comoving frame. Based on the assumption of radiative balance in the fluid rest frame, all ipackets, macro atoms and kpackets have to be converted insitu back to rpackets. For kpackets this is done by sampling the rates at which different physical processes cool the electron gas. All treated cooling rates are listed in Sect. 2.2.3. For macro atoms and ipackets, the situation is more complicated due to the possibility of internal transitions. In both cases we sample the internal energy flow rates until a radiative deexcitaton process is selected or a collisional deactivation to a kpacket occurs (see Lucy 2002, 2003). The needed energy flow rates are calculated with rate coefficients evaluated as described in Sect. 2.2.
2.1.3. Reconstruction of radiation field quantities
For our detailed treatment of ionization and thermal structure (see Sect. 2.2), estimates for the radiative boundfree rates and radiative heating rates are needed. We use volumebased estimators (Lucy 1999b, 2003) to reconstruct the relevant quantities from the trajectories of the packet ensemble. In this approach, the timeaveraged contributions of all trajectory segments, on which the process can in principle occur, are taken into account. Thus, to obtain an estimate for the photoionization rate coefficient γ_{i, j, k} for level i, j, k, we sum over all path segments ds for which the comoving frame (CMF) frequency ν′ of the packet is larger than the threshold for photoionization ν_{i, j, k}
Here, V is the volume of the respective grid cell and is the CMF packet energy. The time interval Δt is a numerical normalization factor that is determined by the energy injection rate at the lower computational boundary. Similarly, the estimator for the stimulated recombination rate coefficient is given by
Here, the Saha factor
enters, which connects the LTE level populations to the ground state population of the next higher ionization stage. The heating rate coefficient for photoionization is
Finally, the heating rate H^{ff} due to inversebremsstrahlung is calculated using
with the freefree opacity χ_{ff}(ν) treated according to Eq. (1). Before concluding our presentation of the reconstruction of radiation field quantities, we stress again that currently the estimators for the boundfree processes γ_{i, j, k}, and are only used for hydrogen.
2.1.4. Relativistic transfer
For photosphericphase SNe II the emergent continuum radiation is created in regions well below the photosphere. In these optically thick regions, the radiation field is essentially isotropic in the fluid rest frame and relativistic frame transformations can significantly modify the energy transport in the ejecta by introducing small anisotropies in the lab frame intensity (see, e.g., Hauschildt et al. 1991; Baron et al. 1996b).
To include relativistic effects in the Monte Carlo simulations, TARDIS uses a mixedframe approach. Radiation–matter interactions are handled in the comoving frame whereas the packet propagation is carried out in the lab frame. Whenever necessary we transform the relevant packet properties between the frames. Compared to Kerzendorf & Sim (2014), we have refined the treatment of relativity by including frame transformations of opacities as well as angle aberration. To transform packet energies and frequencies between observer and comoving frame, we use the full Doppler factor instead of a first order approximation. Expressions for the relevant transformation laws can be found in Mihalas & Mihalas (1984) or, specifically for spherical geometries, in Castor (1972). To be consistent with the adopted frame transformations, the distance to the next possible line interaction is now calculated based on the full Dopplershift formula. As a result, the commondirection frequency surfaces, that is, the surfaces that emit line radiation at the same frequency in the observer frame, are no longer planes perpendicular to the line of sight but have a more complicated geometry as described by the relativistic Sobolev theory of Jeffery (1995).
2.2. Plasma state
The original implementation of TARDIS only features approximate excitation and ionization treatments and a very simplified calculation of the thermal structure. We have considerably refined the determination of the plasma state to adapt the code to SNe II. In particular, we have implemented a full NLTE treatment of excitation and ionization for hydrogen and we employ a thermal balance calculation to infer the temperature structure of the envelope.
The calculation of the plasma state involves a simultaneous determination of the excitation and ionization state of the material as well as the thermal structure. To reduce the complexity of this nonlinear problem, we decouple the solution of the excitation and ionization balance as follows: given an initial guess for the kinetic temperature T_{e} and the electron density n_{e}, we calculate level population fractions as outlined in Sect. 2.2.1. Based on the obtained excitation state, we solve for the ionization balance as described in Sect. 2.2.2. Finally, we compute heating and cooling rates (see Sect. 2.2.3), which are needed for the determination of the thermal structure. An outer iteration loop establishes consistency between excitation and ionization and adjusts the temperature such that thermal balance is enforced (see Sect. 2.2.4).
2.2.1. Excitation
TARDIS offers excitation treatments with different levels of sophistication. Level population fractions f_{i, j, k} = n_{i, j, k}/N_{j, k} can be calculated from the Boltzmann excitation equation, a nebular modification thereof (see Abbott & Lucy 1985) or from the steadystate equations of statistical equilibrium.
For the NLTE excitation calculation, electron number densities have to be specified. In this case, the statistical equilibrium equations for the total system decouple and can be solved for each atomic species individually. In the NLTE treatment of Kerzendorf & Sim (2014), only bound–bound interactions and collisional excitation and deexcitation rates were included. We extend the scheme by including radiative and collisional boundfree processes to obtain a more complete description of hydrogen excitation. The necessary photoionization and recombination rate coefficients γ_{i} and α_{i} are reconstructed from the MC simulation by volume based estimators (see Sect. 2.1.3). The collisional ionization and recombination rate coefficients q_{iκ} and q_{κi} are evaluated according to the approximate formula by Seaton (1962). With these processes included, the rate equation for level i of ion j, k is given by
Here, r_{mi} and r_{im} denote the total rate coefficients at which radiative and collisional transitions between level i and m populate and depopulate level i. In the Sobolev approximation, the rate coefficient for deexcitation from an upper level u to a lower level l is given by
and the excitation rate coefficient is
Here, is the mean intensity at the blue wing of the line, β_{lu} is the Sobolev escape probability (see, e.g., Sect. 4.3.1 of Lucy 2002) and A_{lu}, B_{lu} and B_{ul} are the Einstein coefficients. Electron impact excitation rate coefficients c_{lu} are taken from the approximate formula of van Regemorter (1962) with deexcitation rates evaluated according to detailed balance. For hydrogen levels with principal quantum numbers up to n = 7 we use collision strengths from the detailed ab initio calculations of Przybilla & Butler (2004).
Despite fixing the electron number densities, the system of rates (Eq. (8)) remains nonlinear due to the dependence of the Sobolev escape probabilities on the level populations. We use a standard root finding algorithm to solve for the level population fractions f_{i, j, k} = n_{i, j, k}/N_{j, k} and the ion population ratio N_{j + 1, k}/N_{j, k}^{1}. The convergence of the outer iteration loop that establishes a consistent plasma state is accelerated considerably by including the ion population ratio in the solution of the excitation state.
2.2.2. Ionization
In our detailed treatment of ionization we use the derived level population fractions n_{i, j, k}/N_{j, k} and an initial guess for the electron density n_{e} to calculate the total ionization rate coefficient
and the total recombination rate coefficient
for relevant pairs of ions (j, k), (j + 1, k). We only do this for hydrogen in this work. For all other ions, the Saha factor Φ_{j, k} = (N_{j + 1, k}n_{e})/N_{j, k} and the electron density n_{e} serve as approximations of the total ionization and recombination rate coefficients. The Saha factor Φ_{j, k} is evaluated according to the Saha equation at the local radiation temperature T_{R} or the nebular ionization formula of Mazzali & Lucy (1993), see Eqs. (2) and (3) of Kerzendorf & Sim 2014). Based on these ionization and recombination rate coefficients, we iteratively solve for the ion and electron number densities, assuming ionization equilibrium.
2.2.3. Thermal balance
To complete the description of the plasma state, we need an estimate for the electron temperature T_{e} in the ejecta. In Kerzendorf & Sim (2014), T_{e} was set to 0.9T_{R} following Mazzali & Lucy (1993). We replace this simplified treatment of the thermal structure by a thermalbalance calculation based on the heating and cooling rates of the gas.
The thermal energy content, and therefore the temperature, is determined by the energy exchange between the kinetic energy of the ejecta, the radiative energy pool and the pool of atomic internal energy. This transfer is mediated by adiabatic cooling, collisional transitions as well as boundfree and free–free interactions. Assuming a steadystate, the rates for heating and cooling of the ejecta by these processes must cancel. Thus the electron temperature T_{e} is fixed by the requirement of thermal balance
Here, H^{bf} and C^{fb} denote the rates of heating and cooling by boundfree interactions, H^{ff} and C^{ff} the respective rates for free–free processes. The contributions from collisional excitation, deexcitation, ionization and recombination are C^{exc}, H^{deexc}, C^{ion} and H^{recomb}. The final term C^{ad} describes adiabatic cooling of the envelope due to expansion work.
Specifically, collisional excitations from lower levels l, j, k to level i, j, k remove energy from the thermal pool with a rate
where ϵ_{i, j, k} and ϵ_{l, j, k} are the respective level energies. Correspondingly, collisional ionizations from bound levels of ion j, k contribute
to the total cooling rate.
In turn, atomic internal energy is transfered to the thermal pool by the inverse processes of collisional recombination and deexcitation at rates
and
Thermal electrons moving in the field of an ion j, k emit radiative energy according to (see Osterbrock 1974)
which depends on the ionic charge j − 1, the number density of the respective ion N_{j, k} as well as T_{e} and n_{e}. In addition, energy is continuously removed from the thermal electron pool by radiative recombinations. In terms of the modified rate coefficient
the cooling rate by recombinations to level i, j, k can be written as
Photoionizations, in turn, heat the medium with a rate
Finally, the electron gas continuously loses thermal kinetic due to the expansion of the ejecta. The rate of energy loss resulting from this adiabatic cooling is given by
where t denotes the time of explosion.
2.2.4. Outer plasma iteration
To obtain a consistent solution for the plasma state, the input electron densities that are used in the calculation of the level population fractions have to agree with the results from the ionization calculation. This is achieved by combining the methods described above with an iterative rootfinding procedure. Apart from establishing consistency between the excitation and ionization state, the outer iteration loop is used to determine the thermal structure from the thermal balance Eq. (13).
2.3. Approximations
As established by Utrobin & Chugai (2005), Dessart & Hillier (2007, 2010), Potashov et al. (2017) timedependent effects in the excitation and ionization balance can play an important role in shaping the spectral energy distribution. This applies in particular to epochs following hydrogen recombination. At these times the inclusion of timedependent terms induces an overionization compared to the steady state solution, which is crucial in reproducing the observed Hα line strengths. In contrast, for epochs preceding hydrogen recombination the influence of time dependence is modest. Since these epochs are most relevant for the application of EPM (see, e.g., Dessart & Hillier 2006; Dessart et al. 2008), we do not consider our neglect of these effects a severe limitation to our approach. In fact, Dessart & Hillier (2007) find negligible differences between the dilution factors from their timedependent calculations and the steadystate results from Dessart & Hillier (2006), Dessart et al. (2008). Only for color temperatures less than 7000 K the correction factors drop systematically below their steadystate counterparts. As an additional approximation in the solution of the statistical equilibrium equations, we assume detailed radiative balance in the Lyman continuum. This prevents MC noise in the estimator for the ground state photoionization rate from hindering convergence and avoids associated fluctuations in the ionization as well as the heating and cooling balance. This approach follows previous studies of SNe II, such as Takeda (1990, 1991) and Duschinger et al. (1995). Since the Lyman continuum is optically thick as long as the outflow ionization is not extremely high, detailed balance is deemed to be a very good approximation under most conditions of interest. We have verified this by a series of test calculations without this assumption but with increased numbers of MC quanta. The spectra resulting from the two approaches show good agreement for the parameter space that has been investigated in this paper. This is consistent with the results of Duschinger et al. (1995) who have reached the same conclusion for pure hydrogen supernova atmospheres with photospheric temperatures up to 15 000 K.
2.4. Iteration cycle
TARDIS alternates between the calculation of the plasma state and MC radiative steps to achieve a selfconsistent state for radiation and matter. Generally, less than twenty of these iterations are needed to achieve convergence to a point where only statistical variations remain (see Fig. 1). The good convergence properties result from the strict enforcement of radiative equilibrium, the explicit treatment of scattering and the direct dependence of the macro atom emissivities on the current estimate of the radiation field through the macro atom activation rates. At the moment, the number of iterations is set by hand at the beginning of the simulation, since no formal convergence criterion is implemented in TARDIS. For the calculations presented in this paper, we have performed 40 iterations in all cases. We have found this to be more than sufficient to guarantee convergence for all used setups.
Fig. 1. Illustration of the convergence properties of plasma and radiation field quantities. We show the fractional changes between successive iterations for the mean intensity of the radiation field J, the electron temperature T_{e} and a representative level population, specifically that of the second excited level of hydrogen n_{3}. In all cases, we include both the changes in each individual shell (gray) as well as their average (blue). The results shown here are taken from our SN1999em model for the 14th of November (see Sect. 4 and Fig. 3). 

Open with DEXTER 
2.5. Atomic data
We use the hydrogen atomic data as described by Sim et al. (2005). This data set is based on a 20 level model atom with each level corresponding to a principal quantum number n. Frequencydependent photoionization cross sections are tabulated for every energy state. The tabulated values range from the threshold ionization frequency up to the point at which the crosssection is only about 0.07% of the value at threshold. This improved hydrogen model atom complements the atomic data already included in TARDIS (see Kerzendorf & Sim 2014), which is compiled from the line lists of Kurucz & Bell (1995) and the Chianti 7.1 data base (Dere et al. 1997; Landi et al. 2012).
2.6. Spectral synthesis
To calculate synthetic spectra, the properties of escaping packets are recorded and later binned. However, the quality of the spectra that can be obtained from the normal MC quanta is severely affected by MC noise. To improve the signaltonoise ratio an additional type of MC packet is used in TARDIS. Whenever a normal packet is launched or performs an interaction in the final spectral synthesis run, these so called “virtual” packets (νpackets) are emitted to estimate the contribution of the event to the emergent spectrum. In practice this amounts to optical depth integrations along a number of randomly selected trajectories through the ejecta. The measured optical depths τ_{trj} are then used to weight the contributions of the virtual packets to the spectrum according to the escape probability along the packet path, exp(−τ_{trj}). This procedure, commonly called “peeling off”, is well established and has found widespread use, in particular in the area of dust MC radiative transfer (see, e.g., YusefZadeh et al. 1984; Wood & Reynolds 1999; Baes et al. 2011; Steinacker et al. 2013; Lee et al. 2017).
Compared to the implementation described in Kerzendorf & Sim (2014), modifications have been necessary to keep the computational effort reasonable for the high optical depths of our SN II atmospheres. Since the number of interactions scales quadratically with optical depth, the number of virtual packets that have to be tracked for each “real” packet quickly becomes prohibitively large as the optical depth is increased. At the same time, the contribution of the additional vpackets to the spectrum is marginal due to the strong attenuation towards the surface. We apply biasing to the virtual packet emission to tackle this issue. Virtual packets are created only with a probability exp(−τ_{e}), where τ_{e} is the electron scattering optical depth. To account for the lower chance of creation, the weight of the spawned packet is increased by the inverse of this probability. Notwithstanding this application of biasing, virtual packets can still accumulate large amounts of optical depth, for example, in line interactions. We use the Russian roulette technique (see, e.g., Carter & Cashwell 1975; Dupree & Fraley 2002) to probabilistically remove these lowweight packets.
2.7. Supernova model
TARDIS allows for the use of complex supernova models based on hydrodynamical explosion simulations and with stratified abundances (see Kerzendorf & Sim 2014, Appendix A). Nevertheless, to facilitate the exploration of the parameter space, we restrict our analysis to simple, highly parameterized models. As in D05, we assume powerlaw density profiles
with density indexes n = −dln ρ/dln r in the range n = 6 − 14. Both hydrodynamic simulations (Chevalier 1976, 1982; Blinnikov et al. 2000) and spectral modeling (see, e.g., Eastman & Kirshner 1989; Schmutz et al. 1990; Baron et al. 2007; Dessart & Hillier 2006; Dessart et al. 2008) have demonstrated that the outer density distribution is well described by such an ansatz with values close to n ∼ 10. The composition of the ejecta is taken to be homogeneous. Heavy elements up to nickel are included in the simulations. Following D05 we use CNOcycle equilibrium values from Prantzos et al. (1986) for the abundances of H, He, C, N, O. The remaining elements are assumed to have solar chemical composition with values taken from Asplund et al. (2009).
3. Expanding photosphere method
3.1. Presentation of the method
The expanding photosphere method (EPM) of Kirshner & Kwan (1974) is based on a simplified model of the supernova as a sharplydefined, sphericallysymmetric, expanding photosphere. The radiation emerging from this photosphere is assumed to be that of a blackbody, diluted by an amount given by the dilution factor ξ_{ν}. This correction factor ξ_{ν} has originally been introduced by Hershkowitz et al. (1986a,b) and Hershkowitz & Wagoner (1987) to correct for the dilution of continuum flux that occurs in a scatteringdominated environment. In practice, the dilution factors account for all deviations of the spectrum from blackbody emission, such as lines or limbdarkening, in a parametrized fashion (see, e.g., E96; D05). For reasons of simplicity, in the application of EPM it is assumed that the dilution factor only depends on the color temperature. The precise form of this dependence may be reconstructed from supernovae whose distance is known from independent means (see Schmidt et al. 1992). However, to determine absolute distances it is necessary to infer the dilution factors from theoretical models as in E96 and D05, and outlined at the end of this section.
Based on the assumptions given above, the specific luminosity of the supernova is given by
where R_{ph} is the photospheric radius and T is the temperature of the blackbody B_{ν}(T). By equating this to the observed dereddened luminosity the angular size of the expanding photosphere
can be inferred from the measured dereddened flux . The temperature T has to be determined from photometry as will be outlined shortly. Finally, to obtain the distance to the supernova
the photospheric radius must be eliminated from the equations. For homologous expansion this can be achieved via the relation
where t_{0} is the time of explosion and ν_{ph} is the photospheric velocity. The expansion velocity ν_{ph} can be inferred from blueshift velocities of lines, from crosscorrelation of the observations with model spectra (see Hamuy et al. 2001) or from tailored radiative transfer calculations (see, e.g., Dessart & Hillier 2006; Dessart et al. 2008). Finally, by measuring the ratio of photospheric angular diameter and velocity
for multiple epochs t, the distance is obtained from the slope of the data points. The time of explosion follows from the intercept with the t axis.
To apply this formalism to observations, we have to recast the relevant equations in terms of photometric magnitudes. For a bandpass with a transmission function the apparent magnitude of the object can be calculated from the observed flux according to
where is the zeropoint. Using this definition, we can rewrite Eq. (24) for the diluteblackbody emission as follows:
Here, we have introduced the broadband dust extinction and the blackbody magnitude
where T_{S} is the color temperature. By minimizing the difference between observed and model magnitudes
for a bandpass combination S, the angular diameter θ and the color temperature T_{S} can be inferred from photometric observations.
To determine dilution factors from a synthetic spectrum, Eq. (32) is rewritten in terms of absolute magnitudes :
In this case, the photospheric radius R_{ph} is known and application of the minimization procedure to the synthetic magnitudes yields the color temperature T_{S} and the dilution factor ξ_{S} for the model. In Sect. 6 we will use this procedure to derive an independent set of dilution factors from our TARDIS simulations.
3.2. Dilution factors
To understand the results of our numerical simulations, a firm grasp of the basic physics behind the dilution factors is essential. One of the most important effects in this context and the original motivation for the introduction of the dilution factor (see Hershkowitz et al. 1986a,b) is the dilution of continuum radiation that occurs in a scatteringdominated environment. If, as in SN II, the scattering opacity greatly exceeds the absorptive opacity, a thermally created photon can travel large optical depths before a true absorption event returns it to the thermal energy pool. As a result, these photons can escape the ejecta without thermalizing and can efficiently carry away thermal energy from deep inside the atmosphere. This allows the intensity of the radiation field to fall below the thermal value (B_{ν}) but to still resemble the spectral energy distribution of a blackbody. From random walk arguments, it can be shown (see, e.g., Mihalas 1978) that the relevant optical depth for this process, usually referred to as the thermalization depth τ_{thm}, scales roughly like
where χ_{ν} denotes the total opacity and χ_{ν, abs} the absorptive component. Under these conditions, the emergent flux resembles that of a blackbody with the temperature at the thermalization depth but is diluted by an amount ξ^{2} ≈ 1/τ_{thm}.
4. Example spectra
In Sect. 2, a detailed description of our efforts to extend TARDIS for the spectral modeling of SNe II has been given. Here, we apply the extended code to calculate synthetic spectra for two epochs of SN1999em, a prototypical event of this class. Our goal is to demonstrate that, with the implemented changes, we are able to reproduce the spectral properties of such normal hydrogenrich supernovae. Since we do not aim to perform a quantitative spectroscopic analysis, we have not extensively finetuned the model to exactly fit the observations but have adopted parameters similar to those used in previous studies by Baron et al. (2004) and Dessart & Hillier (2006).
As in Dessart & Hillier (2006), we adopt a powerlaw density profile with index n = 10 and a CNOenhanced composition with an otherwise solar metallicity for both epochs (see Sect. 2.7 for details). Our first model is for the 9th of November, corresponding to around two weeks after explosion. At this point, the hydrogen envelope is still fully ionized but the envelope has already cooled sufficiently for appreciable line blanketing by metals to develop. Apart from the very weak He I 5875 Å feature, helium lines have already disappeared from the spectrum. Since the temperature is still too high for the Ca infrared triplet to form, the spectrum redwards of Hα remains featureless. Our TARDIS model nicely reproduces these characteristics as demonstrated by the comparison to the observations taken by Hamuy et al. (2001) in Fig. 2. The observed spectrum has been dereddened according to a color excess of E(B − V)=0.08, which is slightly less than the value of E(B − V)=0.1 chosen in previous studies by Baron et al. (2004) and Dessart & Hillier (2006). We have blueshifted the observations by 770 km s^{−1} (see Leonard et al. 2002) to correct for the peculiar velocity of the host galaxy. The only major shortcoming of our model is that it underproduces the strength of the Fe II lines at ∼4550 Å and ∼5140 Å. Since this epoch coincides with the recombination from Fe III to Fe II, the predicted strengths of these features are, however, very sensitive to small changes in the parameters and to the adopted ionization treatment.
Fig. 2. TARDIS spectral model (blue) for the observations of SN1999em (black) taken by Hamuy et al. (2001) on the 9th of November. We have smoothed the Monte Carlo spectrum using a Savitzky–Golay filter (Savitzky & Golay 1964). The observational data has been taken from the WISeREP archive (Yaron & GalYam 2012) and has been dereddened according to the Cardelli et al. (1989) law with a color excess of E(B − V)=0.08 (http://www.weizmann.ac.il/astrophysics/wiserep/). To account for the peculiar velocity of the host galaxy, the observations have been blueshifted by 770 km s^{−1} (see Leonard et al. 2002). Finally, the synthetic spectrum has been scaled to match the observed dereddened flux f_{λ}. The main telluric features are marked with circled crosses. 

Open with DEXTER 
The second epoch we are modeling corresponds to an intermediate stage in the photosphericphase evolution of the supernova. On the 14th of November, roughly 3 weeks after explosion, hydrogen recombination has set in and the spectrum shows very prominent Hα emission. Further cooling of the envelope has significantly strengthened the effect of line blanketing compared to the previous epoch. Redwards of Hα the continuum is no longer featureless, since the temperature has dropped sufficiently for the Ca infrared triplet to appear. Figure 3 shows a comparison of the spectrum taken by Hamuy et al. (2001) to our TARDIS model. We have corrected the observations for reddening and peculiar velocities in the same fashion as for the first epoch. Overall, our synthetic spectrum reproduces the measured SED quite well. The two prominent Ca features, Ca H&K and the infrared triplet, are matched well in both strength and shape. However, our model slightly overestimates the Hα emission, whereas the width of the absorption trough is underestimated. As mentioned by Dessart & Hillier (2006), who found similar problems, the latter might be related to blending with Fe II and Si II lines.
Fig. 3. Same as Fig. 2 but for the observations of SN1999em on the 14th of November. 

Open with DEXTER 
5. Model grid
Having established the capabilities of the extended TARDIS version to produce accurate synthetic spectra for SNe II, we use the code to calculate an independent set of dilution factors. To this end, we have set up a grid of 343 models. These have been constructed to cover the interesting physical parameter space of inner boundary temperatures T_{inner} from 9500 to 24 000 K, photospheric densities ρ_{ph} from 7 × 10^{−15} g cm^{−3} to 8 × 10^{−14} g cm^{−3}, powerlaw density indexes n from 6 to 14 and photospheric velocities ν_{ph} from 3000 km s^{−1} to 14 000 km s^{−1}. For these parameters the models span a range of effective temperatures T_{eff} from 4900 K to 12 000 K. We take the photospheric properties ρ_{ph} and ν_{ph} to refer to the position at which the electron scattering optical depth is τ = 2/3. In practice, the models are set up by adopting ρ_{0} = ρ_{*} and r_{0} = ν_{*}t for the density profile (Eq. (23)), where ρ_{*} and ν_{*} are specific values selected from the desired range of photospheric density and velocity. To ensure that the photosphere of the model will lie at the appropriate depth, t (time since explosion) is estimated using
where m_{H} is the mass of a hydrogen atom, μ_{e} is the mean molecular weight per electron and σ_{T} is the Thomson cross section. In making this estimate, it is assumed that μ_{e} is wellapproximated by μ_{e} = 1.52, as appropriate for a composition of ionized hydrogen and singlyionized helium. The full calculation is then carried out and the true values of ρ_{ph} and ν_{ph} are extracted from the simulation. In general, the true values of the photospheric parameters (ρ_{ph}, ν_{ph}) are very close to the originally selected reference parameters (ρ_{*}, ν_{*}) from which the model was generated. Nevertheless, we always refer to each model by the derived (simulation) values of ρ_{ph} and ν_{ph}. We note that the true inner boundary of our computational domain lies considerably deeper than the (approximate) photosphere, typically at τ ∼ 27.
The setup of the model grid is done in the form of a latin hypercube design (Stein 1987). In this approach, each parameter range is subdivided into N equal intervals, where N is the number of models and one random parameter value is selected from each subinterval. This guarantees that, in contrast to a conventional cartesian grid, N distinct values exist for each parameter. This is in particular beneficial if the quantities of interest are only weakly sensitive to a subset of parameters^{2}. Finally, to illustrate the properties of our set of models, a projection of the grid in the T_{eff} – ρ_{ph} plane is shown in Fig. 4. As can be seen, the desired parameter space is for the most part uniformly covered. Deviations from the uniform spacing arise from the use of Eq. (35) to map between photospheric quantities and the computational grid. This process becomes less reliable as soon as a strong recombination front develops. The main motivation for using a mostly uniform grid is that it allows us to study the differential influence of model parameters, such as the photospheric density, on the dilution factors. In this context, correlations between the input parameters have to be avoided as far as possible. However, since quantities such as photospheric temperature and density are certainly not completely independent in nature, this also means that the grid includes models that are not representative of normal SNe II. One example would be an object with a very high expansion velocity but an extremely low temperature.
Fig. 4. Illustration of the model grid from Sect. 5. The plot shows the effective temperature T_{eff} and the photospheric density ρ_{ph} for all atmosphere models. At low temperatures the actual photospheric density can exceed the targeted upper limit of 8 × 10^{−14} g cm^{−3}, due the development of a strong recombination wave, which complicates the mapping between photospheric properties and the computational grid. 

Open with DEXTER 
6. Results
6.1. Overview
From synthetic photometry of our model spectra, we can derive color temperatures T_{S} and dilution factors ξ_{S} according to Eq. (33). To facilitate the comparison to the results of E96 and D05, we focus our analysis on the bandpass combinations S={B,V}, {B,V,I}, {V,I} and {J,H,K} with filter functions taken from Bessell & Brett (1988) and Bessell (1990). Examples of the dilute blackbody models constructed in the synthetic EPM analysis are shown in Fig. 5. The results of our analysis are summarized in Fig. 6, which displays dilution factors ξ_{S} and color temperatures T_{S} for all our models, as well as comparison values from E96 and D05. The color temperatures constitute the most important parameters in the study of the dilution factors, since they account for most of the variance in the correction factors and can be easily inferred from observations. Variations of the remaining parameters such as photospheric density or velocity are in most cases of secondary importance and are responsible for the observed dispersion around the colortemperature trend. In Fig. 6 we find good agreement with D05, in particular at low to medium color temperatures. For {B,V,I} and {V,I} the results match well for temperatures below 12 500 K and for {J,H,K} for temperatures below 7000 K. In {B,V} the dilution factors are similar to D05 over the entire temperature range. For higher color temperatures in {B,V,I}, {V,I}, and {J,H,K} our models tend to be systematically more dilute than D05 with values closer to those published by E96. As will be discussed in Sect. 7.2, part of this discrepancy can be attributed to differences in the adopted photospheric densities. We note that for all bandpass combinations the intrinsic scatter of our dilution factors is slightly larger than for the set of models by D05. This was to be expected, since we have constructed our grid of models in such a way that at all temperatures the whole range of remaining parameters is covered (see Sect. 5). In contrast, the dilution factors of E96 show a much smaller dispersion, since only a narrow part of the parameter space is explored in their study. Before concluding our discussion, we stress that this scatter does not correspond to the diversity of real objects but only reflects our ignorance about the parameter space occupied by SNe II. Finally, following E96 and D05 we present third order polynomial fits to the color temperature dependence of our dilution factors ξ_{S} = ∑_{i}a_{i}(10^{4} K/T_{S})^{i} in Table 1.
Fig. 5. Comparison of the dilute blackbody models from {V,I} (red) and {B,V,I} synthetic photometry (green) to the original TARDIS spectrum (blue). The spectrum is our SN1999em model for the 9th of November (see Fig. 2). The parameters of the blackbody fits are T_{VI} = 9500 K, T_{BVI} = 10 200 K and ξ_{VI} = 0.49, ξ_{BVI} = 0.45. We have overplotted the transmission curves for the B, V and I filters from Bessell (1990). The Monte Carlo spectrum has been smoothed using a Savitzky–Golay filter (Savitzky & Golay 1964). 

Open with DEXTER 
Fig. 6. Dilution factors ξ_{S} as a function of color temperature T_{S} for filter combinations S={B,V}, {B,V,I}, {V,I} and {J,H,K}. We use a common yaxis scale for all bandpass combinations S to highlight the differences in the scaling behavior of the dilution factors. For comparison purposes, we include the models of D05 as blue crosses. Polynomial fits to the dilution factors are shown for all sets of models (dashed: TARDIS, dashed dotted: E96, dotted: D05). For {J,H,K} the D05 curve is not included due to a misprint in the tabulated fit coefficients in the original paper. 

Open with DEXTER 
Coefficients a_{i} of polynomial fits ξ_{S} = ∑_{i}a_{i}(10^{4}K/T_{S})^{i} to the models for bandpass combinations S = {B,V}, {B,V,I}, {V,I}, and {J,H,K}.
6.2. Influence of atmospheric properties
In the previous section, the dependence of the dilution factors on temperature has been discussed. Here, we will focus on the effects of density. First, we analyze the variation of the dilution factors with photospheric density. Secondly, we investigate the influence of the steepness of the density profile. Finally, to conclude the discussion of the effect of model parameters, we will assess the robustness of the dilution factor fit curves against changes in the metallicity.
6.2.1. Influence of the photospheric density
Variations in the photospheric density account for most of the dispersion of the models around the general color temperature dependence as illustrated in Fig. 7. Since the density cannot be easily constrained from observational data, it is essential to understand its influence on the dilution factors to quantify the associated uncertainties on the mean and the variance of the tabulated fit curves. Figure 8 shows the density dependence of our dilution factors and a comparison to the results of D05. In general the dilution factors tend to increase with photospheric density, with the strength of the scaling varying between bandpass combinations. This behavior can be understood by remembering that the amount of continuum flux dilution depends on the ratio of continuum to scattering opacity (see Sect. 3.2). Since the main contribution to the scattering opacity comes from Thomson scattering, it is proportional to the electron density. Thermalization processes, on the other hand, roughly scale with the square of the electron density. As such, we expect the ratio of the two to increase with density, yielding a smaller amount of flux dilution at high densities. To study this behavior in a more quantitative way, we adopt the same ansatz as E96 for the dilution factors:
Fig. 7. Dilution factors ξ_{BVI} as a function of color temperature T_{BVI}. To illustrate the density dependence of our TARDIS models, the logarithm of the photospheric density log_{10}ρ_{ph} is colorcoded. For comparison purposes we include the polynomial fits to the dilution factors of E96 (dasheddotted) and D05 (dashed). 

Open with DEXTER 
Fig. 8. Variation of the dilution factors ξ_{S} with photospheric density ρ_{ph} for filter combinations S = {B,V}, {B,V,I}, {V,I} and {J,H,K}. For comparison purposes the models of D05 are shown as blue crosses. 

Open with DEXTER 
Here, z(T_{S}) denotes a polynomial of the same form as those used in Table 1. A least squares fit to our set of models yields the density scaling indexes γ listed in Table 2. Overall, the inferred density dependence of our dilution factors is moderate with similar magnitudes for all passbands. The scaling indexes γ are systematically a bit larger than those published in E96 with the largest difference in the infrared (see Table 2). However, as illustrated by Fig. 8 the density dependence in the infrared is not well described by a single powerlaw for the entire density range. In Sect. 7.2 we will use the inferred powerlaw scalings to assess whether differences in the assumed photospheric densities play an import role in understanding the discrepancies between the published sets of dilution factors.
Coefficients γ of polynomial fits to the density dependence of the dilution factors for bandpass combinations S = {B,V}, {B,V,I}, {V,I}, and {J,H,K}.
6.2.2. Influence of the density structure
For a powerlaw atmosphere the ratio of photospheric density and density index n is approximately given by
For a given outflow ionization, time of explosion and expansion velocity, an increase in the density index n results in higher photospheric densities ρ_{ph} and therefore less flux dilution. However, if the density structure is treated as an independent parameter, the dilution factors do not show a strong dependence on n as illustrated by Fig. 9. For a possible explanation, we refer the reader to E96, who have found the same behavior and have proposed a physical motivation in their Sect. 3.3.
Fig. 9. Dilution factors ξ for the bandpass combination {V,I} as a function of the density indexes n = −dlnρ/dlnr of the powerlaw model atmospheres. 

Open with DEXTER 
6.2.3. Influence of metallicity
Line blanketing by metals, in particular iron group elements, plays an important role in shaping the spectral energy distribution of SNe II and the resulting influence of metallicity on the emergent spectrum has been discussed in detail in the literature (see, e.g., Dessart & Hillier 2005b). However, neither E96 nor D05 discuss in depth how this effect translates into changes in the dilution factors. To investigate the sensitivity of the ξ − T fit curves to changes in metallicity, we have rerun a random subset of 68 models of our solar metallicity grid (Z = 1) with a lower metallicity of Z = 0.2. The resulting changes in the color temperatures and dilution factors are shown in Fig. 10 for the {B,V,I} bandpass combination. As expected, at high temperatures (T_{BVI} ⪆ 10 000 K) the influence of metallicity on the model properties is negligible, since the degree of ionization is too high for significant line blanketing to develop in the optical and nearUV^{3}. For moderate and low temperatures large changes in the color temperature up to thousands of degrees are observed. However, the associated changes in the dilution factors are approximately aligned with the general scaling of ξ with T. Compared to the intrinsic scatter of the models, the induced changes in the functional behavior of ξ with T are of secondary importance. Thus, in the investigated regime, ranging from solar to distinctly subsolar, the tabulated fit curves are robust against modifications of the metallicity.
Fig. 10. Change of dilution factors ξ_{BVI} with metallicity Z. The subset of the original solar metallicty models for which corresponding calculations with a lower metallicity have been performed are shown in red. The metalpoor models are depicted in blue. Arrows illustrate the changes in model properties induced by the change to the subsolar metallicity of Z = 0.2. To facilitate the comparison to the general ξ − T trend, all models of the original grid are included in gray. 

Open with DEXTER 
7. Comparison to previous studies
7.1. Radiative transfer
To put our results into context, we review differences in the radiative transfer modeling between the CMFGEN code used by D05, the EDDINGTON code used by E96 and TARDIS and discuss possible effects on the dilution factors. The main differences lie in the ionization treatment of metal species, the handling of line opacity and the inclusion of relativistic transfer effects.
For the calculations presented by D05, only the effect of the Doppler shift on the frequency of the radiation field is taken into account. E96 follow a different approach based on the premise that radiationfield time dependence can be included in a quasistatic treatment by enforcing a constant luminosity in the comoving frame. In this case, the timedependent comovingframe transport equation reduces to a much simpler expression that differs from D05 only by an additional term βI_{ν}/r, where I_{ν} is the specific intensity of the radiation field. This term is formally identical to the part of the full transport equation that describes the redshift of photons in the scattering process and thus the adiabatic loss of radiation energy. However, the sign is changed and the magnitude decreased by a factor of three. Both approaches neglect the socalled advection term that arises from the frame transformation of angles (see, e.g., Pistinner & Shaviv 1994, for a discussion). This term is generally deemed to be more important than the aberration term (see, e.g., Baron et al. 1996b). Taking into account the additional reduction of the magnitude of the aberration term in E96, we conclude that the differences in handling the relativistic terms between E96 and D05 are small in comparison to our relativistic treatment (Sect. 2.1.4), which corresponds to a full solution of the quasistatic relativistic transport problem. Since we can achieve good agreement with D05 despite this difference, we consider it unlikely that relativistic effects play an important role in explaining the systematic offset between E96 and D05.
Another possible source for discrepancies, which has been discussed previously in the literature (see D05), is the treatment of line interactions. Here, the differences start with the handling of the opacity. Both CMFGEN (D05) and TARDIS treat the contributions of all lines to the opacity individually, in a consistent manner. In contrast, EDDINGTON E96 adopts the more convenient but approximate expansion opacity formalism of Eastman & Pinto (1993) that combines all line opacity in a wavelength bin. For the opacity calculation, the expansion opacity formalism in E96, as well as the method used in TARDIS, rely on the Sobolev approximation (Sobolev 1957), whereas CMFGEN adopts the comovingframe method. For microturbulent velocities of less than 100 km s^{−1}, as adopted in D05, the Sobolev method is of similar accuracy as the comovingframe method in describing the formation of hydrogen lines in SNe II (see Duschinger et al. 1995). In regions where line overlap is possible, in particular in the metal line forest in the blue, the Sobolev approximation may, however, be less accurate than the comovingframe method.
With respect to line interactions, the final difference between the codes concerns the redistribution of the absorbed radiation. Only CMFGEN computes a full NLTE source function for all included species. In E96 line interactions are treated in detail only for a few selected elements, in most cases only hydrogen. For the remaining species, resonance scattering is assumed and effects such as fluorescence or collisional deexcitations are neglected. TARDIS strikes a balance between the approximate treatment of E96 and the full NLTE calculation of D05. In principle, our implementation of the macro atom scheme of Lucy (2002, 2003) also provides a full NLTE description of the redistribution process. However, in our simulations only radiative and collisional bound–bound transitions are included for species other than hydrogen. Despite this simplification, the approximate NLTE emissivities from the macro atom provide a full treatment of fluorescence. It is also worth pointing out that the predicted emissivities are largely insensitive to errors in the excited states population and therefore to our use of approximate excitation treatments. This has been demonstrated by Lucy (2002) and may be understood by considering that, in the context of the macro atom machinery, the most relevant level number densities are those of the ground state and lowlying metastable levels. Radiative excitations from these states account for most of the activations of the macro atom and their populations are likely to be close to LTE with respect to the ground state. In contrast, the level number densities of excited states, which will be less accurately estimated, are not as important in setting the rate of macro atom activations and enter in the emissivity only through minor modifications of the internal redistribution probabilities for stimulated emission. Thus we argue that the macro atom approach captures most of the essential physics of a full NLTE treatment as opposed to the resonance scattering approximation used in E96. As such, it constitutes a promising source of systematic discrepancies between E96 on the one hand and D05 and TARDIS on the other hand.
To conclude our discussion of the major differences in the numerical treatments, we compare the different methods used for calculating the ionization state of metal species. An accurate solution to the ionization balance is essential in modeling the line blanketing, which shapes the spectral energy distribution in the blue. Due to the use of superlevels, CMFGEN (D05) is able to consistently, though approximately, treat all species in NLTE. In contrast, E96 calculated the ionization using NLTE only for a few selected species and only for a subset of their atmospheric models. For the remaining models and species, LTE at the electron temperature is assumed. Similarly, TARDIS relies on simplified prescriptions for the calculation of the ionization balance of metals. For the results presented in this paper, we have used the nebular ionization approximation of Mazzali & Lucy (1993). In principle, this method should provide a more accurate description of the ionization balance in a diluted, radiation dominated environment than the assumption of LTE. However, neither assumption can fully replace a detailed photoionization calculation. Still, our spectral models for SN1999em (see Sect. 4) reproduce the observed line blanketing well. This instills confidence that, at least for the early and intermediate stage evolution, the nebular ionization treatment adequately captures the essential physics.
Ultimately, it is extremely difficult to assess the extent to which, if at all, individual numerical differences contribute to the systematic discrepancy between the sets of dilution factors. Based on qualitative arguments, we have deemed it unlikely that the handling of relativistic terms plays an important role in this context. We have identified the use of the very simple resonance scattering approximation by E96 as one of the main distinguishing features from both our and D05’s numerical approaches. As such, it can be regarded as a promising possible contributory factor to the systematic differences. However, these interpretations are speculative and should be taken with a grain of salt.
7.2. Effect of model grid assumptions
In the previous section we have discussed how differences in the radiative transfer calculations can affect the dilution factors. In the context of the discrepancy between E96 and D05 most of the discussion in the literature (see, e.g., Dessart & Hillier 2005a; Jones et al. 2009) has revolved around these issues. However, another possibly important (albeit banal) source of systematic differences is the choice of model grid properties. To demonstrate this, we have modified the plot depicting the temperature dependence of our dilution factors for the {B,V,I} bandpass combination to include the colorcoded photospheric density for each model (see Fig. 7). From this, it is obvious that the inferred fit curves can easily be moved upwards or downwards by preferentially sampling either the high density or the low density regions of the parameter space. Since the exact distribution and correlation of parameters such as density, temperature and velocity are not known for the population of SNe II, there exists a certain amount of freedom in the setup of the model grid.
To quantitatively illustrate the role such effects may have, we have investigated the influence of density in particular. For a comparative study we need to consider families of models for which color temperature, densities and dilution factors are available – accordingly, we make use of our TARDIS models, the E96 models and the models presented for the tailored EPM analysis of SN1999em, SN2005cs and SN2006bp (Dessart & Hillier 2006; Dessart et al. 2008). This set of 38 models covers the relevant range of color temperatures and generally follows the fit curves published in D05^{4}.
Before we compare densities, we approximately correct for changes of the electron densities between the set of models due to differences in composition (specifically, we rescale the densities from E96’s models with the estimated ratio of the mean molecular weights per electron). The mean rescaled densities ⟨ρ_{ph}⟩ are shown in Fig. 11 as a function of {B,V,I} color temperature. Overall, the densities used in this paper, and in Dessart & Hillier (2006) and Dessart et al. (2008) tend to be larger than those of E96 with maximum differences of a factor of a few. The conspicuous jump in density for the E96 models between 8000 and 9000 K stems from two exponential atmospheres (e12.2, e12.3). To check whether this density mismatch might alleviate some of the tension between the dilution factors by E96, and those of Dessart & Hillier (2005a, 2006) and Dessart et al. (2008), we rescale the dilution factors ξ_{S} using the simple powerlaw relation from Sect. 6.2.1. For this purpose we do not use the mean densities ⟨ρ_{ph}⟩ from Fig. 12 but the appropriate average . Figure 12a illustrates the effect of the rescaling on the discrepancy between the two sets of models. Applying the density correction reduces the maximum difference from roughly 40% to 20%, but fails to remove the systematic offset completely. As can be seen in Fig. 12b, the procedure is more successful for our set of dilution factors. After rescaling, only a maximum difference of around 8% remains between our calculations and those of Dessart & Hillier (2006) and Dessart et al. (2008). We stress that due to the simplifying assumptions we have made the results above are only qualitative in nature. Nevertheless, our discussion demonstrates that differences in the setup of the model grid, for example, different choices for the photospheric densities, can introduce systematic uncertainties on the 10% level in the dilution factors. This most likely explains part of the discrepancy between the results of E96, and those of Dessart & Hillier (2005a, 2006) and Dessart et al. (2008). To eliminate this additional error source, approaches are needed that strongly constrain the relevant parameters through observational data. One possibility would be to base the dilution factor fit curves on the tailored EPM analyses of a representative set of SNe II.
Fig. 11. Comparison of the mean photospheric density ⟨ρ_{ph}⟩ at a given {B,V,I} color temperature for the models from E96 (red dashed), the tailored EPM analyses of Dessart & Hillier (2006) and Dessart et al. (2008; D06/D08, blue solid), and this paper (green dashdotted). The plotted densities have been rescaled slightly to account for differences in the composition as outlined in Sect. 7.2. 

Open with DEXTER 
Fig. 12. Comparison of the discrepancy between the dilution factor set in the {B,V,I} bandpass combination of Dessart & Hillier (2006) and Dessart et al. (2008) ⟨ξ_{D06}⟩, and E96 ⟨ξ_{E96}⟩ (panel a), and with respect to the results of this paper ⟨ξ_{Tardis}⟩ (panel b). The findings before (blue solid) and after (red dashed) the application of a density correction factor are shown. The details of the procedure are described in Sect. 7.2. 

Open with DEXTER 
8. Conclusions
In this work, we present an extension of the Monte Carlo radiative transfer code TARDIS to the spectral synthesis of SNe II. The key feature of our numerical approach is an updated radiation–matter interaction scheme, which provides a full treatment of bound–bound, boundfree, free–free and collisional processes based on the macro atom scheme of Lucy (2002, 2003). The second major improvement concerns the calculation of the plasma state. The code now contains a selfconsistent determination of the thermal structure from the heating and cooling balance as well as a full NLTE calculation of the ionization and excitation state for hydrogen. Other changes include an improved handling of relativistic effects, an adaption of the spectral synthesis calculation for high optical depths and a different initialization of the plasma state. We demonstrate the capabilities of the extended code by modeling two different epochs of the prototypical SN II SN1999em. For both epochs good agreement with the observed spectra is achieved, instilling confidence that TARDIS is wellsuited for quantitative spectroscopic analysis of photosphericphase, hydrogenrich supernovae.
In line with our goal to use TARDIS for measuring distances, our final application is the calculation of an independent set of EPM dilution factors. In this context, a longstanding issue has been the systematic discrepancy of around 20% between the results of E96 and D05, which translates into an uncertainty of the EPM distance of the same magnitude. To address this problem, we have performed radiative transfer calculations for a set of 343 TARDIS models, which span a wide range of temperatures, densities and expansion velocities. Despite using significantly different numerical techniques, the dilution factors extracted from these calculations show good agreement with those published by D05. This result helps remove some of the tension between the available sets of distance correction factors. It is still somewhat unclear which differences in the numerical approach make the models of E96 systematically more dilute than ours and D05’s. Based on our calculations, we can plausibly rule out only one of the previously suggested explanations, namely differences in the treatment of relativistic effects.
Our other focus lay on investigating the parameter dependences of the dilution factors. Similar to E96 and D05, we identify density as one of the most important parameters in setting the magnitudes of the dilution factors. Our powerlaw fits to the density dependence yield similar scaling behaviors as for the calculations by E96. As in E96, we do not find a strong effect of the steepness of the density profile. In addition, we have demonstrated that changing the metallicity from solar to decidedly subsolar (Z = 0.2) only induces minor modifications in the relationship between color temperature and dilution factors.
Finally, we have investigated differences in the setup of the model grid as an additional source of systematic errors. In our discussion, we have demonstrated that part of the discrepancy between E96 and D05 can plausibly be tracked back to differences in the assumed photospheric densities. This result highlights the need to base tabulated dilution factors on approaches that constrain the model parameters and their correlations more strongly through observational data. One way to achieve this would be to apply the tailored EPM (Dessart & Hillier 2006; Dessart et al. 2008) to a representative set of SNe II.
In this paper we have established TARDIS as a new independent numerical tool for modeling SNe II and have demonstrated its capability to calculate accurate dilution factors. As a next step we plan to apply the code to measure absolute distances using the tailored EPM (Dessart & Hillier 2006; Dessart et al. 2008) or SEAM (Baron et al. 1995, 1996a, 2004, 2007). As a consequence of the inclusion of a much more detailed treatment of the radiative transfer process, the typical runtime of the TARDIS spectral synthesis procedure has increased from minutes needed in the original implementation by Kerzendorf & Sim (2014) to hours. However, in light of the ubiquity of machine learning techniques and the continuous increase in computational resources, this increase in computational complexity is of minor concern and it will, for the first time, be feasible to perform the spectral fitting process in an automated manner. In combination with sampling techniques the parameter space can be explored in a systematic manner and uncertainties in the estimated parameters can be obtained. This will put strong constraints on the accuracy of absolute distance measurements of SNe II and will help to assess their promise as tools for cosmology.
Specifically, we use a modified version of the Powell hybrid method as implemented in SCIPY (Jones et al. 2001).
Acknowledgments
CV thanks Andreas Floers, Stefan Lietzau and Markus Kromer for stimulating discussions during various stages of this project. The authors gratefully acknowledge Stefan Taubenberger for sharing his enormous expertise in the field of supernovae and for being a key member of the supernova cosmology project that motivated this work. CV would also like to thank Michael Klauser for providing help and guidance during the start of this project. The authors thank the anonymous reviewer for valuable comments. This work has been supported by the Transregional Collaborative Research Center TRR33 “The Dark Universe” of the Deutsche Forschungsgemeinschaft and by the Cluster of Excellence “Origin and Structure of the Universe” at Munich Technical University. SAS acknowledges support from STFC through grant, ST/P000312/1. Data analysis and visualization was performed using MATPLOTLIB (Hunter 2007), NUMPY (Oliphant 2006) and SCIPY (Jones et al. 2001).
References
 Abbott, D. C., & Lucy, L. B. 1985, ApJ, 288, 679 [NASA ADS] [CrossRef] [Google Scholar]
 Anderson, J. P., Gutiérrez, C. P., Dessart, L., et al. 2016, A&A, 589, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Baes, M., Verstappen, J., De Looze, I., et al. 2011, ApJS, 196, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Barna, B., Szalai, T., Kromer, M., et al. 2017, MNRAS, 471, 4865 [NASA ADS] [CrossRef] [Google Scholar]
 Baron, E., Hauschildt, P. H., Branch, D., et al. 1995, ApJ, 441, 170 [NASA ADS] [CrossRef] [Google Scholar]
 Baron, E., Hauschildt, P. H., Branch, D., Kirshner, R. P., & Filippenko, A. V. 1996a, MNRAS, 279, 799 [NASA ADS] [CrossRef] [Google Scholar]
 Baron, E., Hauschildt, P. H., & Mezzacappa, A. 1996b, MNRAS, 278, 763 [NASA ADS] [CrossRef] [Google Scholar]
 Baron, E., Nugent, P. E., Branch, D., & Hauschildt, P. H. 2004, ApJ, 616, L91 [NASA ADS] [CrossRef] [Google Scholar]
 Baron, E., Branch, D., & Hauschildt, P. H. 2007, ApJ, 662, 1148 [NASA ADS] [CrossRef] [Google Scholar]
 Bessell, M. S. 1990, PASP, 102, 1181 [NASA ADS] [CrossRef] [Google Scholar]
 Bessell, M. S., & Brett, J. M. 1988, PASP, 100, 1134 [NASA ADS] [CrossRef] [Google Scholar]
 Blinnikov, S., Lundqvist, P., Bartunov, O., Nomoto, K., & Iwamoto, K. 2000, ApJ, 532, 1132 [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]
 Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245 [NASA ADS] [CrossRef] [Google Scholar]
 Carter, L., & Cashwell, E. 1975, Particletransport Simulation with the Monte Carlo Method, Tech. rep. (Los Alamos: Los Alamos National Laboratory) [CrossRef] [Google Scholar]
 Castor, J. I. 1972, ApJ, 178, 779 [NASA ADS] [CrossRef] [Google Scholar]
 Chevalier, R. A. 1976, ApJ, 207, 872 [NASA ADS] [CrossRef] [Google Scholar]
 Chevalier, R. A. 1982, ApJ, 259, 302 [NASA ADS] [CrossRef] [Google Scholar]
 D’Andrea, C. B., Sako, M., Dilday, B., et al. 2010, ApJ, 708, 661 [NASA ADS] [CrossRef] [Google Scholar]
 de Jaeger, T., GonzálezGaitán, S., Anderson, J. P., et al. 2015, ApJ, 815, 121 [NASA ADS] [CrossRef] [Google Scholar]
 de Jaeger, T., Galbany, L., Filippenko, A. V., et al. 2017, MNRAS, 472, 4233 [NASA ADS] [CrossRef] [Google Scholar]
 Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., & Young, P. R. 1997, A&AS, 125, 149 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dessart, L., & Hillier, D. J. 2005a, A&A, 439, 671 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dessart, L., & Hillier, D. J. 2005b, A&A, 437, 667 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dessart, L., & Hillier, D. J. 2006, A&A, 447, 691 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dessart, L., & Hillier, D. J. 2007, MNRAS, 383, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Dessart, L., & Hillier, D. J. 2010, MNRAS, 405, 2141 [NASA ADS] [Google Scholar]
 Dessart, L., & Hillier, D. J. 2011, MNRAS, 410, 1739 [NASA ADS] [Google Scholar]
 Dessart, L., Blondin, S., Brown, P. J., et al. 2008, ApJ, 675, 644 [NASA ADS] [CrossRef] [Google Scholar]
 Dessart, L., Hillier, D. J., Waldman, R., & Livne, E. 2013, MNRAS, 433, 1745 [NASA ADS] [CrossRef] [Google Scholar]
 Dupree, S. A., & Fraley, S. K. 2002, A Monte Carlo Primer (New York: Kluwer Academic/Plenum Publishers) [CrossRef] [Google Scholar]
 Duschinger, M., Puls, J., Branch, D., Hoeflich, P., & Gabler, A. 1995, A&A, 297, 802 [NASA ADS] [Google Scholar]
 Eastman, R. G., & Kirshner, R. P. 1989, ApJ, 347, 771 [NASA ADS] [CrossRef] [Google Scholar]
 Eastman, R. G., & Pinto, P. A. 1993, ApJ, 412, 731 [NASA ADS] [CrossRef] [Google Scholar]
 Eastman, R. G., Schmidt, B. P., & Kirshner, R. 1996, ApJ, 466, 911 [NASA ADS] [CrossRef] [Google Scholar]
 Gall, E. E. E., Kotak, R., Leibundgut, B., et al. 2016, A&A, 592, A129 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gall, E. E. E., Kotak, R., Leibundgut, B., et al. 2018, A&A, 611, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gutiérrez, C. P., Anderson, J. P., Hamuy, M., et al. 2017a, ApJ, 850, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Gutiérrez, C. P., Anderson, J. P., Hamuy, M., et al. 2017b, ApJ, 850, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Hamuy, M., & Pinto, P. A. 2002, ApJ, 566, L63 [NASA ADS] [CrossRef] [Google Scholar]
 Hamuy, M., Pinto, P. A., Maza, J., et al. 2001, ApJ, 558, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Hauschildt, P., Best, M., & Wehrse, R. 1991, A&A, 247, L21 [NASA ADS] [Google Scholar]
 Hershkowitz, S., & Wagoner, R. V. 1987, ApJ, 322, 967 [NASA ADS] [CrossRef] [Google Scholar]
 Hershkowitz, S., Linder, E., & Wagoner, R. V. 1986a, ApJ, 303, 800 [NASA ADS] [CrossRef] [Google Scholar]
 Hershkowitz, S., Linder, E., & Wagoner, R. V. 1986b, ApJ, 301, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Hicken, M., Friedman, A. S., Blondin, S., et al. 2017, ApJS, 233, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres: An Introduction to Astrophysical Nonequilibrium Quantitative Spectroscopic Analysis (Princeton, UK: Princeton University Press), 930 [Google Scholar]
 Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
 Jeffery, D. J. 1995, ApJ, 440, 810 [NASA ADS] [CrossRef] [Google Scholar]
 Jerkstrand, A., Fransson, C., Maguire, K., et al. 2012, A&A, 546, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python, http://www.scipy.org/ [Google Scholar]
 Jones, M. I., Hamuy, M., Lira, P., et al. 2009, ApJ, 696, 1176 [NASA ADS] [CrossRef] [Google Scholar]
 Kerzendorf, W. E., & Sim, S. A. 2014, MNRAS, 440, 387 [NASA ADS] [CrossRef] [Google Scholar]
 Kirshner, R. P., & Kwan, J. 1974, ApJ, 193, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Kurucz, R. L., & Bell, B. 1995, Atomic Line List (Cambridge: Smithsonian Astrophysical Observatory) [Google Scholar]
 Landi, E., Del Zanna, G., Young, P. R., Dere, K. P., & Mason, H. E. 2012, ApJ, 744, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, G. K. H., Wood, K., DobbsDixon, I., Rice, A., & Helling, C. 2017, A&A, 601, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lentz, E. J., Baron, E., Branch, D., & Hauschildt, P. H. 2001, ApJ, 557, 266 [NASA ADS] [CrossRef] [Google Scholar]
 Leonard, D. C., Filippenko, A. V., Gates, E. L., et al. 2002, PASP, 114, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Lucy, L. B. 1999a, A&A, 288, 282 [NASA ADS] [Google Scholar]
 Lucy, L. B. 1999b, A&A, 345, 211 [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]
 Magee, M. R., Kotak, R., Sim, S. A., et al. 2016, A&A, 589, A89 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Magee, M. R., Kotak, R., Sim, S. A., et al. 2017, A&A, 601, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mazzali, P. A., & Lucy, L. B. 1993, A&A, 279, 447 [NASA ADS] [Google Scholar]
 Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. (San Francisco: W. H. Freeman and Co.) [Google Scholar]
 Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics (New York: Oxford University Press), 731 [Google Scholar]
 Mitchell, R. C., Baron, E., Branch, D., et al. 2002, ApJ, 574, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Oliphant, T. E. 2006, Guide to NumPy (Provo, UT: Brigham Young University) [Google Scholar]
 Osterbrock, D. E. 1974, Astrophysics of Gaseous Nebulae (San Francisco: W. H. Freeman and Co.) [Google Scholar]
 Pistinner, S., & Shaviv, G. 1994, ApJ, 432, 302 [NASA ADS] [CrossRef] [Google Scholar]
 Potashov, M. S., Blinnikov, S. I., & Utrobin, V. P. 2017, Astron. Lett., 43, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Poznanski, D., Butler, N., Filippenko, A. V., et al. 2009, ApJ, 694, 1067 [NASA ADS] [CrossRef] [Google Scholar]
 Prantzos, N., Doom, C., de Loore, C., & Arnould, M. 1986, ApJ, 304, 695 [NASA ADS] [CrossRef] [Google Scholar]
 Przybilla, N., & Butler, K. 2004, ApJ, 609, 1181 [NASA ADS] [CrossRef] [Google Scholar]
 Quimby, R. M., Wheeler, J. C., Hoflich, P., et al. 2007, ApJ, 666, 1093 [NASA ADS] [CrossRef] [Google Scholar]
 Rodriguez, O., Clocchiatti, A., & Hamuy, M. 2014, AJ, 148, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Savitzky, A., & Golay, M. J. E. 1964, Anal. Chem., 36, 1627 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, B. P., Kirshner, R. P., & Eastman, R. G. 1992, ApJ, 395, 366 [NASA ADS] [CrossRef] [Google Scholar]
 Schmidt, B. P., Kirshner, R. P., Eastman, R. G., et al. 1994, ApJ, 432, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Schmutz, W., Abbott, D. C., Russell, R. S., Hamann, W.R., & Wessolowski, U. 1990, ApJ, 355, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Seaton, M. J. 1962, in Atomic and Molecular Processes, ed. D. R. Bates (New York: Academic Press), 375 [Google Scholar]
 Sim, S. A., Drew, J. E., & Long, K. S. 2005, MNRAS, 363, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Sobolev, V. V. 1957, Sov. Astron., 1, 678 [NASA ADS] [Google Scholar]
 Stein, M. 1987, Technometrics, 29, 143 [CrossRef] [Google Scholar]
 Steinacker, J., Baes, M., & Gordon, K. D. 2013, ARA&A, 51, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Takats, K., & Vinko, J. 2006, MNRAS, 372, 1735 [NASA ADS] [CrossRef] [Google Scholar]
 Takeda, Y. 1990, A&A, 234, 343 [NASA ADS] [Google Scholar]
 Takeda, Y. 1991, A&A, 245, 182 [NASA ADS] [Google Scholar]
 Utrobin, V. P., & Chugai, N. N. 2005, A&A, 441, 271 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Regemorter, H. 1962, ApJ, 136, 906 [NASA ADS] [CrossRef] [Google Scholar]
 Wood, K., & Reynolds, R. J. 1999, ApJ, 525, 799 [NASA ADS] [CrossRef] [Google Scholar]
 Yaron, O., & GalYam, A. 2012, PASP, 124, 668 [NASA ADS] [CrossRef] [Google Scholar]
 YusefZadeh, F., Morris, M., & White, R. L. 1984, ApJ, 278, 186 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Coefficients a_{i} of polynomial fits ξ_{S} = ∑_{i}a_{i}(10^{4}K/T_{S})^{i} to the models for bandpass combinations S = {B,V}, {B,V,I}, {V,I}, and {J,H,K}.
Coefficients γ of polynomial fits to the density dependence of the dilution factors for bandpass combinations S = {B,V}, {B,V,I}, {V,I}, and {J,H,K}.
All Figures
Fig. 1. Illustration of the convergence properties of plasma and radiation field quantities. We show the fractional changes between successive iterations for the mean intensity of the radiation field J, the electron temperature T_{e} and a representative level population, specifically that of the second excited level of hydrogen n_{3}. In all cases, we include both the changes in each individual shell (gray) as well as their average (blue). The results shown here are taken from our SN1999em model for the 14th of November (see Sect. 4 and Fig. 3). 

Open with DEXTER  
In the text 
Fig. 2. TARDIS spectral model (blue) for the observations of SN1999em (black) taken by Hamuy et al. (2001) on the 9th of November. We have smoothed the Monte Carlo spectrum using a Savitzky–Golay filter (Savitzky & Golay 1964). The observational data has been taken from the WISeREP archive (Yaron & GalYam 2012) and has been dereddened according to the Cardelli et al. (1989) law with a color excess of E(B − V)=0.08 (http://www.weizmann.ac.il/astrophysics/wiserep/). To account for the peculiar velocity of the host galaxy, the observations have been blueshifted by 770 km s^{−1} (see Leonard et al. 2002). Finally, the synthetic spectrum has been scaled to match the observed dereddened flux f_{λ}. The main telluric features are marked with circled crosses. 

Open with DEXTER  
In the text 
Fig. 3. Same as Fig. 2 but for the observations of SN1999em on the 14th of November. 

Open with DEXTER  
In the text 
Fig. 4. Illustration of the model grid from Sect. 5. The plot shows the effective temperature T_{eff} and the photospheric density ρ_{ph} for all atmosphere models. At low temperatures the actual photospheric density can exceed the targeted upper limit of 8 × 10^{−14} g cm^{−3}, due the development of a strong recombination wave, which complicates the mapping between photospheric properties and the computational grid. 

Open with DEXTER  
In the text 
Fig. 5. Comparison of the dilute blackbody models from {V,I} (red) and {B,V,I} synthetic photometry (green) to the original TARDIS spectrum (blue). The spectrum is our SN1999em model for the 9th of November (see Fig. 2). The parameters of the blackbody fits are T_{VI} = 9500 K, T_{BVI} = 10 200 K and ξ_{VI} = 0.49, ξ_{BVI} = 0.45. We have overplotted the transmission curves for the B, V and I filters from Bessell (1990). The Monte Carlo spectrum has been smoothed using a Savitzky–Golay filter (Savitzky & Golay 1964). 

Open with DEXTER  
In the text 
Fig. 6. Dilution factors ξ_{S} as a function of color temperature T_{S} for filter combinations S={B,V}, {B,V,I}, {V,I} and {J,H,K}. We use a common yaxis scale for all bandpass combinations S to highlight the differences in the scaling behavior of the dilution factors. For comparison purposes, we include the models of D05 as blue crosses. Polynomial fits to the dilution factors are shown for all sets of models (dashed: TARDIS, dashed dotted: E96, dotted: D05). For {J,H,K} the D05 curve is not included due to a misprint in the tabulated fit coefficients in the original paper. 

Open with DEXTER  
In the text 
Fig. 7. Dilution factors ξ_{BVI} as a function of color temperature T_{BVI}. To illustrate the density dependence of our TARDIS models, the logarithm of the photospheric density log_{10}ρ_{ph} is colorcoded. For comparison purposes we include the polynomial fits to the dilution factors of E96 (dasheddotted) and D05 (dashed). 

Open with DEXTER  
In the text 
Fig. 8. Variation of the dilution factors ξ_{S} with photospheric density ρ_{ph} for filter combinations S = {B,V}, {B,V,I}, {V,I} and {J,H,K}. For comparison purposes the models of D05 are shown as blue crosses. 

Open with DEXTER  
In the text 
Fig. 9. Dilution factors ξ for the bandpass combination {V,I} as a function of the density indexes n = −dlnρ/dlnr of the powerlaw model atmospheres. 

Open with DEXTER  
In the text 
Fig. 10. Change of dilution factors ξ_{BVI} with metallicity Z. The subset of the original solar metallicty models for which corresponding calculations with a lower metallicity have been performed are shown in red. The metalpoor models are depicted in blue. Arrows illustrate the changes in model properties induced by the change to the subsolar metallicity of Z = 0.2. To facilitate the comparison to the general ξ − T trend, all models of the original grid are included in gray. 

Open with DEXTER  
In the text 
Fig. 11. Comparison of the mean photospheric density ⟨ρ_{ph}⟩ at a given {B,V,I} color temperature for the models from E96 (red dashed), the tailored EPM analyses of Dessart & Hillier (2006) and Dessart et al. (2008; D06/D08, blue solid), and this paper (green dashdotted). The plotted densities have been rescaled slightly to account for differences in the composition as outlined in Sect. 7.2. 

Open with DEXTER  
In the text 
Fig. 12. Comparison of the discrepancy between the dilution factor set in the {B,V,I} bandpass combination of Dessart & Hillier (2006) and Dessart et al. (2008) ⟨ξ_{D06}⟩, and E96 ⟨ξ_{E96}⟩ (panel a), and with respect to the results of this paper ⟨ξ_{Tardis}⟩ (panel b). The findings before (blue solid) and after (red dashed) the application of a density correction factor are shown. The details of the procedure are described in Sect. 7.2. 

Open with DEXTER  
In the text 