Shock cooling emission from explosions of red super-giants: II. An analytic model of deviations from blackbody emission

Light emission in the first hours and days following core-collapse supernovae (SNe) is dominated by the escape of photons from the expanding shock heated envelope. In a preceding paper, Paper I, we provided a simple analytic description of the time dependent luminosity, $L$, and color temperature, $T_{\rm col}$, valid up to H recombination ($T\approx0.7$ eV), for explosions of red supergiants with convective polytropic envelopes without significant circum-stellar medium (CSM). The analytic description was calibrated against"gray"(frequency-independent) photon diffusion numeric calculations. Here we present the results of a large set of 1D multi-group (frequency-dependent) calculations, for a wide range of progenitor parameters (mass, radius, core/envelope mass ratios, metalicity) and explosion energies, using opacity tables that we constructed (and made publicly available), including the contributions of bound-bound and bound-free transitions. We provide an analytic description of the small, $\simeq10\%$ deviations of the spectrum from blackbody at low frequencies, $h\nu<3T_{\rm col}$, and an improved (over Paper I) description of `line dampening' for $h\nu>3T_{\rm col}$. We show that the effects of deviations from initial polytropic density distribution are small, and so are the effects of `expansion opacity' and deviations from LTE ionization and excitation (within our model assumptions). A recent study of a large set of type II SN observations finds that our model accounts well for the early multi-band data of more than 50\% of observed SNe (the others are likely affected by thick CSM), enabling the inference of progenitor properties, explosion velocity, and relative extinction.


INTRODUCTION
In core collapse supernovae (SNe) explosions, a radiation mediated shock (RMS) traverses outwards through the stellar progenitor, heating and expelling material as it passes.If no significant circumstellar material (CSM) is present around the star, arrival of the shock at the surface produces a hard UV/X-ray ∼ 10 45 erg s −1 'shock-breakout' emission, lasting from tens of minutes to an hour.The breakout pulse is then followed in the coming hours and days by thermal UV/optical 'shock-cooling' emission, caused by diffusion of photons out of the shock-heated stellar ejecta.Typical luminosities and temperatures during shock-cooling are of the order 10 42 − 10 44 erg s −1 , and 1 − 10 eV.As the photons diffuse out, deeper parts of the ejecta are gradually exposed over time (see Waxman & Katz 2016;Levinson & Nakar 2020, for reviews).
In order to constrain the properties of the progenitor star, it is helpful to have high cadence multi-band observations in the first hours of shock-cooling (see Morag et al. 2023, for a detailed discussion of the theoretical considerations and observational status).Among these measurements, ultraviolet observations are especially important, as they are closer to the thermal emission peak, and can ★ E-mail: jonathan.morag@weizmann.ac.il be used to determine the emission temperature and the UV extinction self-consistently (Rabinak & Waxman 2011;Sapir & Waxman 2017;Sagiv et al. 2014;Rubin & Gal-Yam 2017).Combined with an accurate theoretical model, these measurements can be used to reproduce the progenitor and explosion parameters, including radius, surface composition, explosion energy per unit mass, and the extinction.Analytic models are especially important for solving the "inverse problem" of inferring system parameters and uncertainties from the observed spectral energy distribution (SED).
Emission during shock-cooling is amenable to modeling in the case of 'envelope breakout' (i.e. in the absence of CSM with signifcant optical depth), largely because the system is near local thermal equilibrium (LTE) at this time.However, catching supernovae within the first hour presents a practical challenge, and few such observations have been achieved (see Morag et al. 2023, for a representative list).Existing and upcoming observatories, such as the Zwicky Transient Factory (ZTF) (Gal-Yam et al. 2011), the upcoming Vera Rubin Observatory (Ivezić et al. 2019), and the expected launch of the wide-field UV space telescope ULTRASAT (Sagiv et al. 2014;Shvartzvald et al. 2023) will greatly increase the quantity and quality of early measurements, enabling a systematic study.
The numeric calculations presented here provide a systematic analysis of the deviation of the spectra of the emitted radiation from blackbody, which enables an improvement of the accuracy of the analytic models.We adopt a 'multigroup' (MG) numeric approach, where radiative transfer is solved under the diffusion approximation including coupling to hydrodynamics.Several frequencydependent codes using different approximations and schemes have been used to address the shock breakout and cooling problem: STELLA (Blinnikov & Bartunov 2011) is a 1 dimensional code that solves for the angle-dependent intensity with a variable Eddington factor, employing a ray-tracing scheme.Their opacity includes free-free, bound-free and atomic lines from Kurucz (1995), and employs the Sobolev approximation, based on Eastman & Pinto (1993).STELLA should contain the required physical components for modeling the SED of shock-cooling emission, and this has been done in several works (Blinnikov et al. 1998(Blinnikov et al. , 2000;;Tominaga et al. 2011Tominaga et al. , 2009;;Förster et al. 2018).Sapir & Halbertal (2014) numerically solved the planar shock breakout problem under the diffusion approximation and incorporating free-free opacity.They also include inelastic Compton scattering, which is important for shock breakouts at high velocities.Other hydrodynamically coupled multigroup codes are presented in the literature (Vaytet et al. 2011;Skartlien 2000;Tominaga et al. 2015), though they have not been used or formulated to solve for shock-cooling emission.
Many other codes focus computational efforts on line modeling and are specialized to describe emission at later phases of the supernova, when the ejecta is freely-coasting.The ARTIS (Kromer & Sim 2009) and SEDONA (Roth & Kasen 2015) codes both solve frequency-dependent, time-dependent radiative transfer using Monte Carlo.ARTIS is primarily used to solve type Ia SN problems, while SEDONA has also been used to calculate core-collapse supernova emission (e.g.Kleiser & Kasen 2012;Tsang et al. 2020;Jacobson-Galán et al. 2022).The latter contains a module for coupling hydrodynamics to radiation in a Lagrangian sense, though to our knowledge this has not been used in the context of the first day of shock-cooling emission.Other multigroup codes solve for radiative transfer using static "snapshots" of the hydrodynamic profiles, under the assumption that the photon diffusion time is short compared to the dynamical time of the ejecta (Baron et al. 2000;Dessart et al. 2015;Baron et al. 2003;Kerzendorf & Sim 2014;Ergon et al. 2018).It is not guaranteed that this approximation is valid at the earliest times, though our results suggest that it may be reasonable close to recombination if one employs realistic density and temperature profiles.
In a preceding paper, Morag et al. (2023, hereafter Paper I), we modeled shock-cooling emission following breakout from a spherical envelope, assuming local thermal equilibrium (LTE).Using earlier analytic results (Sapir et al. 2011;Katz et al. 2012;Sapir et al. 2013;Rabinak & Waxman 2011;Sapir & Waxman 2017, hereafter SKW I, KSW II, SKW III, RW11, SW17), we derived an analytic description of the bolometric luminosity and color temperature col , that provides a good approximation (to 10% in and 5% in col ) of the results of our 'gray' diffusion simulations for the emission over ∼ 1 hour to ∼ 1 week after shock breakout for a wide range of explosion energies and progenitor parameters (radii, core and envelope masses, metallicities).The analytic expression is advantageous relative to those of previous analytic works (SKW I; KSW II; SKW III; RW11; Nakar & Sari 2010;Shussman et al. 2016;Piro et al. 2017) thanks to its calibration using a large set of numeric results incorporating a realistic opacity with free-free, bound-free, and bound-bound components, which have an important effect on emission.We also used a set of MG diffusion simulations to show that the spectrum is well described by a Planck spectrum with an effective color temperature, except at UV frequencies where the flux may be significantly suppressed due to line absorption.
In this work we relax the assumption that the photons are at LTE and account for frequency-dependent emission, absorption and transport both numerically and analytically.We produce hydrodynamically coupled 1-dimensional multigroup simulations for a wide range of explosion energies and progenitor parameters that span the parameter range of Paper I -explosion energies in the range = 10 50 − 10 52 erg, progenitor masses and radii = 2 − 40 ⊙ , = 3 × 10 12 − 2 × 10 14 cm, envelope to core mass ratios env / c = 10 − 0.3, and metallicities = 0.1 ⊙ − ⊙ .Our simulations are based on the work of Sapir & Halbertal (2014).Their code allows inelastic Compton scattering, which we do not incorporate in the bulk of the simulations, due to its negligible effect during shock-cooling.We include a frequency-dependent opacity , for which we use both the publicly available TOPS code (Colgan et al. 2016) and one that we developed ourselves and is now available to the community Morag (2023).Relative to TOPS, our code is largely open source, and can produce opacity tables at arbitrary density, temperature, and frequency resolution.It is based on experimentally verified atomic line lists (Kurucz 1995).We note that our MG simulations cannot be used to make predictions for the exact spectrum including lines, as they have finite frequency resolution, and do not include a treatment of "expansion opacity" effects and of deviations of excitation and ionization states from LTE (see § 7.1).However, they are useful for predicting the coarser-grained, line averaged SED in shock-cooling following envelope breakout.
A recent study of a large set of type II SN observations (Irani et al. 2022) finds that our model accounts well for the early multi-band data of 50% of observed SNe, corresponding to 70% of the intrinsic SN distribution (the others are likely affected by thick CSM).This agreement enables the inference of progenitor properties, explosion velocity, and relative extinction by using the formulae provided in this paper.
This paper is organized as follows.In § 2 we define our notation and summarize the analytic results of SKW I, KSW II, RW11, SW17 and Paper I, that we use in this paper.In § 3 we describe our numeric calculations and present convergence and code verification tests.In § 4 we describe our opacity tables, as well as convergence and verification tests for these.In § 5 we derive an analytic description of the deviations of shock cooling emission spectra from blackbody, and in § 6 we compare the analytic description to simulation results.In sections § 7.1 and § 7.2, we address the sensitivity of our results to "expansion opacity" corrections and to deviations of plasma ionization and excitation from those of LTE.A comparison to earlier work, including to STELLA radiation transport calculations for several non-polytropic profiles obtained using MESA stellar evolution calculations (Blinnikov et al. 1998;Tominaga et al. 2011;Kozyreva et al. 2020), is given in § 8.The agreement of our results with these earlier calculations provides additional support to our code's validity, and to the conclusion of the detailed analysis of SW17, who demonstrated that the shock cooling emission is not sensitive to deviations of the density profile from a polytropic one.Our results are summarized and discussed in § 9.

RMS BREAKOUT AND SHOCK-COOLING EMISSION-SUMMARY OF EARLIER ANALYTIC RESULTS USED IN THIS PAPER
RMS breakout and shock-cooling is extensively discussed in the literature and briefly summarized in Paper I. At radii close to the stellar radius ( ≡ ( − )/ ≪ 1), the initial density of a polytropic envelope approaches a power-law, 0 = ¯ . (1) Here ¯ ≡ /(4 3 /3) is the average pre-explosion density of the ejecta (exculding the mass of a possible remnant), and = 3/2 for convective RSG envelopes.is a numerical factor, of order unity for convective envelopes, that depends on the inner envelope structure (Matzner & McKee 1999;Calzavara & Matzner 2004;Sapir & Waxman 2017).The predicted breakout and cooling emission are nearly independent of .
As the shock approaches the edge of the star, it accelerates down the steep density profile and the flow approaches the self-similar solutions of Gandel 'Man & Frank-Kamenetskii (1956); Sakurai (1960).The shock velocity diverges in this regime as with 1 = 0.19, and with v s * a constant defined by eq.(2).Based on numerical calculations, Matzner & McKee (1999) find where is the mass of the ejecta, is the energy deposited in the ejecta, and v * is its characteristic expansion velocity.This approximation holds to better than 10% for env / c < 1/3, and overestimates v s * by approx.20% for env / c = 0.1 (see figure 7 of SW17).
Breakout takes place when the scattering optical depth of the plasma layer lying ahead of the shock equals es = /v sh (Ohyama 1963) 1 .We denote the shock velocity v sh and pre-shock envelope density 0 at this point in terms of breakout parameters bo and v bo , respectively.We may rewrite eqs.( 1) and (2) as The location at which breakout "occurs", i.e.where = /v sh , is given by bo where is the opacity.For RSGs, bo and v bo are approximately related to the progenitor parameters and explosion energy by (Waxman & Katz 2016) Here, = 10 13 13 cm, = 0.34 0.34 cm 2 g −1 , v * = v * ,8.5 10 8.5 cms −1 , and = 1 0 ⊙ .The duration over which the breakout pulse bo is emitted from the star is approximately given by the shock crossing time of the breakout layer bo , bo where bo = 10 −9 bo,−9 g cm −3 and v bo = 10 9 v bo,9 cm s −1 .The observed pulse duration may be longer than this intrinsic duration due to light travel time effects, which spread the pulse over / .bo is given as a function of progenitor parameters and explosion energy as bo = 0.02 0.90 13 ( 0 v s * ,8.5 0.34 where v s * = v s * ,8.5 10 8.5 .
For later use, we provide here the density profile during the spherical phase, given by eq. ( 9) of RW11.We recast this equation in terms of and using RW11 eqs.(3), ( 4) and ( 8 In Paper I, we described shock cooling emission by interpolating between the exact planar phase solution (SKW I; KSW II; SKW III) valid at early times (hours), and the later approximate spherical phase solution (RW11; SW17).The combined bolometric luminosity and emission (color) temperature col are given by col / col,br = min 0.97 ˜ −1/3 , ˜ −0.45 .
Assuming a blackbody spectral distribution, the emitted luminosity is then given by, Here { , , } = {0.9, 2, 0.5} (note slight difference from SW17), and ˜ ≡ / br .tr is roughly the time at which the photons will be able to diffuse out of the envelope in dynamical time.The br (break) subscript marks the values at the transition between the planar and spherical phase.
The former is the time past which we showed light-travel time effects to be unimportant.The latter is the time at which the photosphere temperature is = 0.7 eV, based on RW11, roughly corresponding to Hydrogen recombination.The transparency time, tr 2 , occurs roughly when the dynamical time matches the diffusion time, given by SW17 For later use we also define the homologous time hom , which is also the early validity of the RW11; SW17 formula: hom = /5 s, * = 0.1 13 / s * ,8.5 days.
During shock-cooling, the luminosity is determined at the diffusion depth, the location from which photons will diffuse outwards in dynamical time (see e.g.Rabinak & Waxman 2011).The color temperature meanwhile, is roughly determined by the temperature at the thermal depth, the last absorption surface for diffusing photons (the radius from which photons diffuse out of the ejecta without further absorption).In this paper we treat these quantities as frequencydependent (compare with the approximate prescription, Paper I, eq.30).Specifically the thermal depth col is defined by where the abs, es and subscripts indicate absorption, (electron) scattering and frequency dependence, respectively.This integral is often approximated in the literature as 3 abs es or 3 abs ( abs + es ).We find that the choice of approximation has a negligible effect on the SED due to the steep ( ) dependence, with the exception of regimes with strong lines, where the observed effect can be tens of percents.

DESCRIPTION OF THE NUMERICAL CODE
We numerically solve the radiation hydrodynamics equations of a spherically symmetric flow.Following Sapir & Halbertal (2014), we allow the photon distribution to deviate from thermal equilibrium with a multi-group treatment, and handle radiative transfer under the diffusion approximation.For the matter equation of state (EOS), we assume an ideal Hydrogen gas in LTE, including ionization as dictated by the Saha equation.The MG emission/absorption and diffusion opacities (effective opacities for each photon energy group), are based on our tables with a solar mix-like composition.This section is structured as follows.The equations of the numerical scheme are given in § 3.1 and § 3.2, and the initial and boundary conditions are described in § 3.3.The validation of the numeric code and the convergence of the calculations are described in § 3.4.

Radiation-hydrodynamics equations
Our equations employ the diffusion approximation and omit terms of order ( / ) 2 for radiative transfer, based on Castor (2007).In Lagrangian coordinates, the velocity v and density evolve in response to the radiation energy density per unit frequency and matter energy density , as follows: 2 Our 0.7 eV and tr correspond roughly to 2 and 4 (or 5 ) of Utrobin (2007).
where is the plasma pressure, * is the freq.dept.diffusion opacity, is the photon flux, and the 0 subscript denotes initial values.Correspondingly, the (Lagrangian, fluid element associated) energy densities evolve according to, The emission / absorption, compression, and diffusion terms are given by emis/abs where the frequency-dependent photon energy flux density is given self-consistently by and the Planck energy density per unit frequency is given by 4 (34) * is the diffusion opacity, given by the Rosseland mean of sum of the Thomson scattering opacity es and the absorption opacity abs, across the frequency bin (of each photon energy group in the MG calculation), , where abs, is given by our high-resolution opacity tables.The absorption opacity of each group is given by abs = (Δ bin ) −1 ∫ abs, .Eq. ( 29) is based on the assumption that the emitting plasma is in LTE (see e.g.Pinto & Eastman 2000, and further discussion in § 9).
In a small fraction of the simulations we also include ineslastic Compton scattering using the Kompaneets equation, as shown below.Inelastic Compton scattering is unimportant during shock-cooling, but can have an important effect on the breakout spectrum, which we include in our tests of the code in § 3.4.
where e is the electron mass and is the matter temperature (in units of energy).The flux is initially started as 0 = −1 es 3 0 .The time derivative term in eq. ( 33) becomes important in optically thin regions.In our problem it does not have an important effect since the luminosity is determined deep in the ejecta.Following Sapir & Halbertal (2014), we rewrite the radiative compression term in eq. ( 30) as We include Hydrogen recombination in the EOS, given by where ( ) is the ionization fraction, is the Rydberg energy, is the atomic number density, = 5/3 is the monotonic adiabatic index.This prescription becomes an ideal gas equation of state when the plasma is fully ionized.is solved for iteratively using the Saha equation, assuming the presence of Hydrogen only.

Numerical scheme
In our numerical scheme, we solve the continuity equations (eqs.24-26) by a standard staggered mesh leap-frog method.Energy evolution in time is solved via operator splitting.The equations are divided into parts, with diffusion (including a flux limiter), radiative processes and compression calculated consecutively as follows.
Frequency dependent diffusion, is solved implicitly using a Newton Raphson (NR) solver.The output, including and , is then fed into the radiative processes (emission, absorption, and scattering), which are solved iteratively using two loops of NR solvers that each solve for energy conservation without updating .In the inner loop, we solve implicitly for , while keeping constant, while in the outer loop we solve implicitly for , including matter compression.The inner loop includes several protections from non-physical results in , including an 'overshoot' protection to prevent from crossing ( ) when attempting to equalize = ( ) (see details in § 3.2.2).The initial guess for involves solving radiative transfer explicitly, and if the NR solver fails to find a solution after 30 iterations, the solver also attempts a solution starting from the original value.The energy evolution step finishes with radiative compression and then is updated again via eq.33.
Finally, where needed (e.g.eqs.26, 29, 33, * , abs, ), we extract the temperature and pressure by solving the equation of state, eq. ( 37) implicitly.The entire set of equations for the evolution of the matter and radiation energies is solved using a predictor-corrector with opacities updated at every iteration prior to the diffusion step.

Time steps
The minimum of the following constraints limits our simulation time step.For grid cells , the usual Courant upper limit is Δ c = c min Δ i / s,i , where Δ is the grid spacing and is the speed of sound.Diffusion also limits the time step according to , where the minimum is taken over cells and bins centered at frequency .The factor d , along with all the similar factors here, is of order unity and smaller than one, and our results are shown to be insensitive to the exact value.Finally, we limit the time step also by limiting the maximal change due to radiative processes of the total energy density of the radiation/plasma, is the frequency-dependent absorption opacity, is the blackbody distribution, and T is the matter temperature.The radiation temperature is defined as in Sapir & Halbertal (2014) as

'Protections' on Diffusion and Radiative Transfer
The very high opacity, reaching ∼ 10 6 cm 2 g −1 at some frequencies, may lead to numeric problems in the application of the implicit solution with finite time steps.At infinitely small time steps, will be kept close to at such large opacity regions.However, using finite time steps, the implicit result for can at times 'overshoot' (or proceed in the wrong direction due to strong dependence of on temperature).In order to avoid impractically short time steps, we add several limiting 'protections' immediately after the inner Newton Raphson solver for .Namely, if the resulting implicit lies outside of the range between of the previous step and , we override the result to the nearest of these values (see also, Mihalas & Mihalas 1999;McClarren et al. 2008;Gasilov et al. 2016).We also limit the change during the emission/absorption time step to Δ abs,nu / < abs , where abs is in the range {0.1, 0.5}.Its value does not affect our results.Though the latter constraint affects the relative rates of physical processes, absorption still proceeds quickly relative to the other processes in this scheme.
We also add a flux limiter to the simulations.The P1 diffusion approximation doesn't require one in principle.However, at certain frequencies, strong absorption and the strong sensitivity of to temperature, can lead to a situation where the photon energy density in a particular frequency group and cell depletes abruptly -for example to match ( )-.The change may occur faster than the time in which diffusion can respond given finite time resolution, and thus may lead to non-physical flow.Namely, either flux would flow in the wrong direction, or | | > cell , where cell is the energy density in the cell from which the flux is exiting, and is the speed of light.In lieu of adding another time constraint, we insert a flux limiter as follows, where is a positive integer constant.The derivative for -defined in between cells-with respect to u, in either one of the adjacent cells, is written as where ≡ / and ≡ / cell .This derivative is used in the Newton Raphson solver during the diffusion step.
We test the flux limiter on the gray diffusion runs with various values of , finding a deviation of less than a percent from the nonflux limited runs when = 8, which we choose in our simulations.

Initial and boundary conditions
Our numeric calculation involves a succession of three simulations, with each simulation starting later in physical time using a snapshot of the hydrodynamic variables as described by its predecessor.Each successive simulation also contains increasing physical complexity; Tests of our multigroup code.Top: an example of a radiation mediated shock (RMS) in a uniform medium, with photon energy density / ℎ , shown at several locations across the shock.This is compared with an analytic Wien profile, defined by the analytic bolometric energy density , and the local photon number density , extracted from simulation.The shock is of velocity =10%, in a pure Hydrogen medium, with initial proton number density of p = 10 15 g cm −3 .Bottom: luminosity / ℎ , at shock breakout (over the course of several shock-crossing times 0 ) with a solar composition and a free-free opacity only.The simulation is compared at breakout time to the table results of SKW III, where we extract breakout velocity, density, and breakout time based on a fit of the bolometric luminosity to SKW I (as is done in Paper I), without any additional fitting.The progenitor radius is = 10 13 cm, core and envelope masses are env = core = 10 ⊙ , and the explosion energy is = 10 52 erg cm −3 .
i.e. first a hydrodynamic-only calculation, then a gray diffusion calculation, and finally a MG simulation.This way, later time stages of interest include all the relevant processes, while allowing the computations to be performed in practical time.All simulations are proceeded through to the latest times for comparison.
Following SW17 and Paper I, we begin with a simplified progenitor structure, comprising of a uniform density core surrounded by a polytropic envelope at hydrostatic equilibrium.We start a hydrodynamic only simulation where we inject a high thermal energy density in the innermost cells of the core and capture the resulting shock using artificial viscosity.Then a radiative diffusion-hydrodynamics gray simulation is started between 24 and 8 shock crossing times prior to shock breakout, with the exact start time have negligible effect in this range.Both of the simulations are identical to the ones in Paper I (and described there in detail), with the important exception that in the gray simulation we increase the initial cell resolution towards the stellar edge 3 .
Then we begin a multi-group diffusion simulation, as described in § 3.2.The simulation is started between 20 shock crossing times prior to breakout and up to 2 R/c times after breakout, with the exact time at which the simulation is started having negligible effect on the later shock-cooling emission.Typical resolutions for each of the respective stages listed above are 4000-8000, 1600-3200, and 200-1600 cells (MG runs that are started after shock breakout typically have 50-200 cells), with 32-256 photon groups in the MG phases.All calculations are continued until at least after the recombination time 0.7eV .
Multigroup simulations that are started prior to breakout have a similar initial grid to the gray diffusion simulations.Namely the grid changes smoothly, with modest resolution in the interior, highest resolution at the starting location of the shock, and steadily decreasing resolution outwards (keeping cell count constant across the RMSsee SW17), before approaching a constant for / bo , and finally increasing resolution at the stellar edge, with Δ ∼ ( − ), down to at least a scattering optical depth of es ∼ 10 −2 .MG simulations that were started after shock breakout have a simpler initial cell grid, spaced logarithmically in es , with the same stellar edge resolution.All MG simulations have photon frequency bins that are constant in time and are spaced logarithmically.
For all simulations we assume a static reflective boundary in the inner surface, and a free boundary at the outer surface that for the diffusion simulations accelerates as b = b / , where the subscript b denotes boundary values.The boundary flux is given by b = edd b , where edd = 0.3 − 0.5 is the Eddington factor.The results are insensitive to the exact value of edd since the flux is determined deep within the plasma, at ∼ / ≫ 1.

Code validation and numerical convergence
In Paper I, we validated our numerical hydrodynamical-only code against the analytic planar stellar breakout solutions of Gandel 'Man & Frank-Kamenetskii (1956) and Sakurai (1960), and our gray diffusion code against the analytic planar "Sakurai-Weaver" Anzats solutions of SKW I. We also reproduced the bolometric breakout flux expected from a planar stellar breakout, as also described in SKW I.In Sapir & Halbertal (2014), an earlier version of the multi-group code that we use here, underwent several test problems involving radiative diffusion, emission/absorption, and inelastic Compton scattering.
We perform here two additional tests of the multigroup code; the problem of a steady planar radiation mediated shock, and the breakout spectrum in a hydrogen dominated stellar envelope, both including inelastic Compton scattering and only free-free absorption opacity.
We calculate the structure of the steady planar radiation mediated shock at two representative velocities, = v/ = 1% and 10%, i.e. spanning breakout velocities in our parameter range.We find for both cases that density , velocity v, and bolometric photon energy density , converge to 3%, 0.5%, and ∼ 1% relative to the analytic result (see Weaver 1976, and SKW I).For the = 1% simulation, the photons are in LTE, and is in excellent agreement with a Planck distribution matching total bolometric luminosity .For = 10%, the photons carrying most of the energy are in Compton equilibrium.While an exact analytical solution does not exist for the spectrum in this case, we find good agreement between the numeric results for the frequency of the peak of the radiation energy density behind the shock and an analytic estimate using a Wien distribution4 , see figure 1.At lower frequencies, the energy density transitions to a thermal distribution due to large free-free opacity, producing a visible deviation from the Comptonized spectrum, as is expected.
Next, we compare our results for envelope breakout with the approximate table values from SKW I; SKW III, which again assume a fully Comptonized Wien spectrum.We find reasonable (10's of %) agreement in peak temperature and luminosity (fig.1), which is somewhat remarkable, since the temperature and luminosity profiles in these tables are a function of only the breakout parameters ( , bo , bo -see § 2).For the comparison we extract these parameters from the simulation without performing additional fitting for the SED.There is again a noticeable deviation in the low energy tail due to thermalization.Sapir & Halbertal (2014) performed MG diffusion calculations of the planar envelope breakout phase, and obtained similar results to ours.They find 10's of % agreement with SKW I; SKW III in peak temperature and luminosity, as well as a thermalized low-frequency tail of similar shape.
In Paper I, we showed convergence of the hydrodynamic and gray simulations.Our MG calculations are also converged with respect to spatial resolution.Doubling the spatial resolution produces at most a few percent change (and often less than 1%) in ≡ / .We also verify that we are converged with respect to the resolution of the outermost cell, finding that varies by less than 1% when the minimum scattering optical depth varies between es ∼ 10 −2 and 10 −3 , in agreement with the conclusions of Tominaga et al. (2011).We note that our frequency resolution is high enough that our SED is insensitive to the number of photon frequency groups, but is coarse relative to the atomic line scale, which we discuss at length in § 7.1.
As described in Paper I, we include in both the gray and MG diffusion simulations, a non-radiating plasma component coupled with artificial viscosity .This addition helps stabilize against numerical instabilities associated with the density inversion that occurs at the outer edge of the ejecta.

OUR COMPOSITE OPACITY TABLE
Calculating the frequency dependent opacity requires employing several approximations and assumptions.Primary challenges involve solving the many-electron Schrödinger equation and estimating microplasma interactions between species.The assumption that all degrees of freedom are in thermal equilibrium is often, though not always, employed.Due to these approximations, there are large uncertainties in the opacity.For example, a factor of 2 discrepancy in the Rosseland mean exists between TOPS and OP (Iglesias & Rogers 1996) in our regime of interest (as shown in Paper I).
We built our own frequency-dependent opacity table, containing free-free, bound-free, and bound-bound components, and assuming local thermal equilibrium.The code that produces these tables is now available to the community on github Morag (2023) and produces tables at arbitrary density, temperature, and frequency resolution (we use Δ / ∼ 10 −5 − 10 −6 in practice).The opacity code does not include "expansion opacity" effects.
In Paper I we used the Rosseland mean opacity from the highresolution tables to provide a formula for blackbody emission.Here we bin the tables in frequency to produce our multigroup opacities.For the absorption term ( abs, in eq.29) we use the average opacity across each bin, while for the diffusion term ( * in eq.32) we take the Rosseland mean across the bin.In this work, we include 10 important atoms up to iron: H, He, C, N, O, Ne, Mg, Si, S, Fe (though the opacity code can handle an arbitrary mixture).We show later that the resulting supernova lightcurves are insensitive to the exact composition in the case of Hydrogen dominated envelopes.
As discussed in Paper I, the frequency-dependent TOPS opacities are verified against experiments at temperatures and densities exceeding tens of eV and 10 −6 g cm −3 respectively (Colgan et al. 2015(Colgan et al. , 2016(Colgan et al. , 2018)).The Kurucz database of atomic transitions was calibrated against measured line frequencies and oscillator strengths (Kurucz 1995), but is incomplete for highly ionized species at high, > 10 eV, temperature.We separately run simulations using both our opacity table and using TOPS.
This section is written as follows.In § 4.1, we describe the construction of the opacity, and in § 4.2 we describe tests and results of the table.

Opacity -Construction
We solve for the ionization and excitation population numbers selfconsistently with the Saha equation, and using electron levels provided by the NIST database 5 .For the free-free components and Hydrogen-like bound-free and bound-bound components, we use the equations provided in Rybicki & Lightman (1979).We verify our free-free calculations against TOPS, finding 7% agreement or better.
For line transitions in Hydrogen-like ions with charge , the decay rate is given by ( ) ∼ 2 = 4 ( = 1), where ( = 1) is the decay rate for the equivalent Hydrogen line, and ∼ 2 is the frequency of the atomic transition.Non Hydrogen-like boundfree absorption is based on the table in Verner et al. (1996)  6 .In general, the results of the table agree w/TOPS up to a factor of 2, with the exception of a pure Fe mix we tested, where we observe an order of magnitude difference.Bound-bound oscillator strengths, degeneracies and lower energy levels are taken from Kurucz CD 237 .
Our line sampling method varies with line width relative to the grid resolution Δ grid .Lines with full width half max (FWHM) that are much thinner than the grid spacing, Δ , are sampled as a single grid point, with magnitude proportional to ∼ (Δ ) −1 to conserve total flux.Broad lines are sampled as a Voight function , including native line width and thermal broadening.oscillator strength is conserved.We find that we are insensitive to the exact choice of cutoff between each of the sampling regimes in the ranges 1/30<FWHM/Δ grid <1/3 and 1<FWHM/Δ grid < 10, respectively.We therefore choose the cutoffs at FWHM/Δ grid = 1/108 and 3.For quicker calculation in the latter two cases, the broad Voight wings are interpolated at 100 times coarser resolution than the frequency grid, and combined at the end of the calculation9 .Finally, ions and electronic states with electron occupation fractions below 10 −14 are removed (we tested that this omission has no effect).
We also include Red-wing continuum suppression (following TOPS - Colgan et al. 2018), and the Hummer-Mihalas factor that suppresses higher electron states (Hummer & Mihalas 1988).Our opacity code allows computation of other ingredients that were not 10 -2 10 0 10 2 10 4 10 0 10 5 Figure 3.An example of our frequency dependent opacity table for a solar mix at a density of = 10 −11 g cm −3 and temperature = 1 eV, compared with TOPS, showing reasonable agreement.The observed discrepancy in the lower-frequency atomic lines is most likely a result of the lower resolution in TOPS, where such lines cannot be seen.Since these lines are very narrow, the effect on the MG simulations is small.
added in this work as they were shown to have negligible influence on the multigroup simulations in our range of parameters, including the Dappen Anderson Mihalas factor (Dappen et al. 1987), and electron collisional broadening.These latter processes are described in the documentation for the opacity.

Tests of the Opacity Calculation and Results
We test our frequency-dependent opacity code against TOPS and OP for the simplest case of a pure Hydrogen mixture, finding good agreement ( 15%) in both Rosseland mean and frequency-dependent opacity (see figs. 2 and 3).As a further sanity check, we also compare our bound-bound opacity for a pure Fe mix with an independent code from Waxman et al. (2018), finding excellent agreement (not shown).We find that our simulations are converged with respect to the underlying opacity table, as tested in more than 5 separate parameter choices spanning progenitor radius and explosion energies for 128 photon groups.We observe < 2% difference in the SED when varying between Δ / ∼ 10 −5 − 10 −6 base grid, and < 1% when changing the number of { ≡ / 3 , } grid points from {16,66} to {30,120}.

DEVIATIONS FROM BLACKBODY -CALIBRATED ANALYTIC MODEL
Here we provide an analytic description of the deviations of the emitted spectrum from blackbody.Our approximation formulae are derived piecewise, based on the strength of the absorption opacity, abs, , relative to the scattering opacity, es .At frequencies where abs, > es , the emitted flux may be approximated as a blackbody with a frequency-dependent thermal depth (surface of last absorption) col, , and corresponding frequency-dependent color temperature col, = ( col, ), given by At frequencies where the absorption opacity is smaller than the scattering opacity, abs, < es , we base our approximation on the flux emitted by a semi-infinite planar slab of temperature in the two-stream approximation (Rybicki & Lightman 1979) where = abs, /( abs, + es ). (44) We therefore suggest as an approximation for the escaping spectral luminosity in this regime (see also, Dessart & Hillier 2005).
We first derive an expression describing the emission at the regime of relatively low absorption opacity, abs, < es , which occurs primarily at intermediate frequencies near and below the Planck peak (e.g.ℎ ∼ 1−3 eV in fig.3).Absorption in this regime is dominated by free-free transition with a small bound-free contribution.Neglecting the bound-free contribution we approximate eq. ( 23), that defines the frequency-dependent thermal depth, as (see also, SWN16) Here we have neglected abs, with respect to es and used where the density is in cgs, temperature is in ergs.We approximate the gaunt factor ff = √ 3 0 (ℎ / ) ∼ 0.717 (ℎ / ) −0.27 , with 0 being the zeroth modified Bessel function of the second kind.Solving eq. ( 46) using the analytic RW11/SW17 spherical phase density profiles (eqs.9 -12) we obtain the radius, temperature and opacity at the thermal depth,  (48) col, = 2.13 0.28 13 ( 0 ) −0.02 v 0.17 ( We find that modifying the expression for col, to col, = +1.29×10The thermal depth values can then be inserted into eq.( 45) with = ff, /( ff, + es ) in order to describe the emission in the low absorption frequency range.
For frequency regions with strong absorption, where abs, < es , we return to eq. ( 42).The thermal depth at these frequencies is located at the outer edge of the ejecta, where the density decreases sharply and the temperature, determined by the free-streaming photons, is nearly uniform.We therefore approximate col, ≈ .( ) and col, ≈ .( ) for these frequencies, and describe the emission as a gray blackbody , eq. ( 15).At low frequencies where the free-free opacity dominates, we find numerically that the luminosity is well approximated by (0.85 col ).Meanwhile, at frequencies near and above the Planck peak, where atomic transitions dominate, we use both the simulations and a separate analytic estimate (see § 7.1) to improve upon the approximate BB (0.74 col ) description of the UV suppression of Paper I, replacing the suppression factor 0.74 with a function of ( , ) lying in the range [0.6, 1].The combined freq-dept formula is thus where = 5, and , is again given by eqs.( 45) and ( 44) with the choice abs,nu → ff, .The 1.2 factor accounts for modest UV excess we observe in our results at the planck peak due to the presence of strong lines.The frequency slope in the Raleigh-Jeans regime is similar, but slightly lower than the blackbody value ∼ 2 .
Eq. ( 56) can be further simplified to be given in terms of only and col , with a minor decrease in the approximation's accuracy, .and col are given by eqs.( 13) and ( 14).

NUMERIC RESULTS
In figs. 4 and 5 we plot shock-cooling numeric results from over a dozen multigroup simulations.Deviations from our gray blackbody formula (eqs.13-15) are small.In the Raleigh Jeans regime, the SED slope is slightly shallower than ∼ 2 .As frequency increases in this regime, passes from 10's of percents above the Planck distribution to 10's of percents or more below it.Near the Planck peak at ℎ ∼ 3.5 eV, the SED can be 10's of % in excess of the prediction of our blackbody formula (50% in extreme cases, further discussed in § 7.1).In the Wien tail, the flux is suppressed due to line absorption.
The SED obtained numerically is also compared to our frequency dependent formula (eq.56) yielding good agreement, with an RMS error of Δ / 20% for ℎ < 3 col (and several 10's of % or more in the Wien tail, reflecting a very small inaccuracy in the radiation temperature).The frequency-dependent formula is generally closer to the simulation results than the gray formula prediction throughout.Both formulas shown in the figure use breakout parameters ( bo , bo , ) that are derived by comparing the breakout bolometric luminosity to SKW I (as was done in Paper I), without additional fitting.
We find in many simulations that the ratio of bolometric luminosity between the multigroup and gray simulations is approximately given by ( es / Ross ), where the Rosseland mean opacity Ross is evaluated at the diffusion depth, where is determined (recall that the gray simulation only include Thomson opacity).This factor is less important during the planar phase, but generally reduces in the multigroup simulation by 10 ′ of % during the spherical phase, primarily due to Wien suppression.
We test the sensitivity of our results to the choice of opacity table using simulations that span the progenitor radius and explosion energy parameter range.When we use TOPS table instead of our own, we find that near recombination time, col in the presence of the TOPS-derived opacity can be lower by up to 10's of % (usually <5% -see example in fig.6).This result is in general agreement with our gray analysis in Paper I, where we concluded that for col < 4 eV, it is preferable to use our opacity table due the presence of lab-confirmed lines.For higher temperatures, early during shock-cooling, we find negligible (few percent) SED difference between the two opacity tables, since most of the observable frequencies (ℎ 10 eV) are in the Rayleigh Jeans regime, and are less affected by the presence of lines.
We also test the sensitivity to composition in our simulations, finding a negligible variation when metallicity varies from = 0.1 solar to solar metallicity (see fig. 6), in agreement with Paper I results.We conclude that the SED in Hydrogen dominated envelopes is insensitive to metallicity.We note, however, that for near zero metalicity, up to a factor of 2 differences may be obtained at UV frequencies.

Expansion opacity and finite frequency resolution
Our numeric calculation method did not include "expansion opacity" effects and assumed that the opacity, , does not vary significantly across the extent Δ of a single spatial resolution element.The large velocity gradients in the flow we are considering, combined with the presence of strong absorption lines, may lead to significant variations of across Δ due to the space dependent Doppler shift.This effect may have a significant impact on the calculation of photon transport.For the diffusion calculation, we use a Rosseland average of the opacity, , = 1/ −1 where the average is over the frequency band -+1 , implying a photon mean free path of = 1/( , ).The Rosseland mean is dominated by the frequencies at which the opacity is low, i.e.where lines are absent, and hence the presence of strong lines does not affect significantly.In the presence of a large velocity gradient, the photon mean free path may become significantly shorter-as the photon propagates through the varying velocity of fluid, the frequencies of the lines shift and the photon will be absorbed when it reaches a position where the shifted frequency of a strong line coincides with the photon's frequency (Friend & Castor 1983;Eastman & Pinto 1993;Castor 2007;Rabinak & Waxman 2011).In other words, the effective Doppler broadening "closes" the low "windows" in frequency space, that allow a low value of and large . The absorption mean free path is given in this case by where Δ is the frequency separation between strong lines, and we have used / = / as appropriate for the homologous expansion phase.A line is considered "strong" if it leads to > 1 when integrated over the photon propagation path taking into account the Doppler shift-for line opacity = ( − ), = ∫ ′ = ( / ) (where ′ is the Doppler shifted opacity).Thus, Δ is the frequency distance between lines with > 1. ( exp can be used to define an effective "expansion opacity", taking into account the finite velocity difference across Δ due to the expansion of the plasma.The ratio between the Rosseland and expansion mean free paths is where / is the (strong) line density and we have used ∝ −10 (as appropriate for the shock cooling expanding plasma) yielding 10 = .In the region where the escaping flux is determined, ≈ /v, we have exp / ≈ 10( /v) 2 (Δ / ) ≈ 10 3.5 (Δ / ).We therefore expect exp / ≫ 1 in this region, and thus only a negligible effect of the "expansion opacity" on the escaping flux.This is demonstrated to be the case in Fig. 7, showing / exp for various photon energies near recombination time as a function of es , the electron scattering optical depth.As demonstrated in the figure, the suppression of the mean free path due to the "expansion opacity" effect at this time may affect the spectrum at photon energies > 5 eV, since exp / < 1 may be obtained at the thermalization depth of such high energy photons.We arrive at the same conclusion from examination of simulations spanning different progenitor radii and explosion energies.
In order to test the possible impact of the expansion opacity on the high energy spectrum, we compare (see fig. 8) the flux obtained in the numeric simulation, which does not include the effects of expansion opacity, with the flux obtained using the analytic approximations of eqs.( 42) and ( 45) with col, (eq.23) calculated with the density and temperature profiles obtained in the simulation.In the latter we use a high resolution, Δ / ∼ 10 −5 opacity table and take into account the Doppler shifts of the lines 10 .The SED is given by where cut determines the scattering optical depth es at the thermal depth below which we neglect the effect of scattering.The results are not sensitive to the choice of cut in the range 0.3-3.We show results for cut = 1.The analytic approximation, with col calculated using a high resolution frequency-grid and including expansion opacity effects, reproduces well the spectrum obtained in the simulations.The modifications of col , and the implied flux modification in the analytic approximation, due to the inclusion of the Doppler shifts of lines using a high resolution frequency grid, yield therefore an estimate  (eqs. 13, 14) is in black, and our modified frequencydependent formula (eq.56) is in dashed colored lines, showing excellent agreement with the simulations.The blue / red / yellow color denotes results at times 3 / , 1.5eV , 0.7eV .In plots with with black triangles, the time 0.7 eV (in yellow) is replaced by the transparency validity time tr / , since the latter occurs prior to recombination.The formula is based on breakout parameters, ( bo , bo , ), extracted from simulation without additional fitting involved. of the magnitude of this effect.The good agreement between the simulations results and the analytic estimate obtained as described above implies that the simulation, with its coarse frequency resolution, reasonably estimates the effect of the expansion opacity with an accuracy of typically 10%.The comparison of the analytic and numeric results described above, reveals several frequency bins at intermediate photon energies, 5-8 eV, in which the flux obtained in the numeric simulation exceeds that which is obtained by the analytic estimate by a factor of a few.This is due to the presence of strong and relatively isolated lines, which lead to a large average absorption opacity abs,i of the photon bin, which in turn leads to a blackbody photon spectrum across the frequency bin with temperature corresponding to the plasma temperature.Since the radiation energy density at intermediate photon energies is below the plasma temperature at radii where the optical depth is small (see § 7.2 and fig.9), the flux obtained in bin is significantly larger than the flux obtained at other frequencies.This result is an artifact of the numeric calculation using finite frequency resolution-the emitted flux would match the blackbody spectrum (of the plasma temperature) only across the (very) narrow line width, hardly affecting the total flux in the finite frequency range of bin .
In order to remove this effect, we have identified and removed a few such isolated lines from the opacity table used in the simulations.The 5 removed lines are listed in table 1.Much of the remaining discrepancy between the simulations and the analytic estimate is likely due to the fact that the flux obtained in the numeric simulation exceeds that of the analytic results due to an overestimate of the effect of lines.

Deviation from LTE Ionization and Excitation
Deviations of plasma excitation and ionization states from LTE ('NLTE' effects) may become pronounced, and have been studied by many groups, at late times, in particular during the nebular phase (Baron et al. 1996;Utrobin & Chugai 2005;Dessart & Hillier 2008;Lisakov et al. 2017).NLTE effects at the very early times that we are considering, prior to H recombination, are less widely studied and are expected to be small (see, e.g.Blinnikov et al. 2000;Baron et al. 2000;Dessart & Hillier 2008)  11 .We briefly explain below why the 11  Takeda (1990Takeda ( , 1991) ) predicts order of magnitude deviations of atomic electron populations from LTE at ph ∼ 1 eV.We note, however, that these calculations are performed for densities, which are orders of magnitude smaller than those relevant for the shock-cooling problem that we are considering.The Takeda 1990 photospheric radius, 10 14 cm, implies a photospheric density of ∼ 10 −13 g/cc (for the steep 1/ 10 density profile expected for the expanding stellar envelope).This is nearly 3 orders of magnitude higher than the density used in the Takeda 1990 calculation.42) and ( 45) with col, (eq.23) calculated with the density and temperature profiles obtained in the simulation, but with a high resolution, Δ / ∼ 10 −5 opacity table and taking into account the Doppler shifts of the lines (The gray approximation, eqs.( 13)-( 15) shown in black).The high-resolution spectrum has been averaged over frequency bins, without affecting the SED.Top (bottom): progenitor radius = 3 × 10 12 cm ( = 10 14 cm), explosion energy = 10 51 erg, and core and envelope masses env = core = 10 ⊙ .
LTE assumption is reasonable for modeling the SED in the first hours and days.
The density of the plasma emitting the radiation during the hours to days of shock breakout and cooling is relatively high, and the radiation is close to thermal equilibrium with the plasma -see fig.9.This situation is quite different from that prevailing later, on weeks time scale, where the density is low and the radiation is far from thermal equilibrium, causing the ionization fraction and excitation level distribution to deviate largely from a Boltzmann distribution.
The relatively large density implies that the time scale for all relevant processes (electron-electron and electron-ion collisions, photoionization and excitation, electron impact excitation), with the exception of electron impact ionization, are much shorter than the dynamical time (see fig. 10).This, combined with the fact that the photon spectrum at energies exceeding the Planck peak (e.g.ℎ 3 eV during recombination), is close to thermal out to very low es , implies that the ionization fraction and the excitation level distribu- tion of the low energy excited states are both close to LTE.The fact that at low optical depth, es 1, the radiation energy density falls below that of LTE at low photon energy, implies that the level distribution of the higher energy excited states may deviate from LTE at the outer edge of the ejecta.However, we note that the distribution of excited energy levels are strongly dominated by UV transitions to and from the highly populated ground state (accounting for roughly half the rate, e.g.Zel'dovich & Raizer 2002), hence deviations from LTE occupation of the higher energy states are expected to be mild.In addition, the effect of lines on the SED in this energy range is small, and hence we do not expect a significant effect due to deviations from LTE.

Shussman et al. (2016) (herafter SWN16
) produced an analytic model including deviations from blackbody due to free-free opacity.They arrive at an SED that disagrees with our simulations by an order of magnitude or more in the infrared and up to a factor of two near the Planck peak.In Kozyreva et al. (2020), the SWN16 formula is compared to simulations produced by the STELLA code.They conclude, similarly, that the SWN16 model's infrared behavior does not describe the SED for ℎ < 3 col , and suggest using a blackbody formula for this range.In fig.11, we compare the SWN16 formula, including the Kozyreva et al. (2020) corrections, with the results of our MG simulations.We find that their analytic model does not reproduce well the simulation results, due the different values obtained in their analytic approximation for and col , which in turn is largely due to their opacity approximation neglecting bound bound transitions (see Paper I).
Next, we compare our numeric results to those obtained using STELLA (Blinnikov & Bartunov 2011) radiation transport calculations (Blinnikov et al. 1998;Tominaga et al. 2011;Kozyreva et al. 2020) for several non-polytropic profiles obtained using MESA stel- lar evolution calculations (Paxton et al. 2018).In each case, we approximate the progenitor density profile used in earlier simulations by a simple polytrope.We show two such examples in fig.12 and summarize the profile parameters we use in table 2. In figs.13 and 14 we then compare the resulting emission from our multigroup simulations to published results, without performing any further fitting.The comparison to the results of Tominaga et al. (2011) is shown in fig.13.We find a good agreement with our results when we choose not to include bound-bound opacity.When the line opacity is turned on, as we normally do for our simulations, there is clear disagreement with regards to the peak shape and temperature, though the Rayleigh-Jeans regime shows excellent agreement.
The result in Tominaga et al. (2011) appears to us as non-physical, since the SED peaks at a wavelength corresponding to ∼ 3 times the color temperature, that is expected based on density and temperature profiles (and on the effective opacity-the scattering optical depth at the location of the corresponding temperature is (2020) correction (dashed gray).At top and bottom, explosion energies of = 10 50 and 10 52 erg, with progenitor radius = 10 13 cm, core and envelope masses of env = core = 10 ⊙ .very high, ≈ 100), and since the SED exhibits a sharp photoionization cutoff, which we do not expect to be observed at this resolution as it should be completely occluded by nearby lines.
The comparison to the results of Kozyreva et al. (2020) is shown in fig.14 (top) during shock breakout.We find an excellent agreement between the results of the different calculations (note that at these temperatures only few atomic transition lines contribute).
Finally, fig.14 (bottom) shows a comparison to the results of Blinnikov et al. (1998).There is reasonable agreement between the results at the latest observed times ( = 3 days after breakout).At early time, there is good agreement at longer wavelengths, > 3000Å, where the free-free emission dominates, but there is a clear discrepancy at higher frequencies.The very high temperature of the emission peak obtained in Blinnikov et al. (1998) at these times, ∼ 70 eV, appears to be non-physical, since similarly to the Tominaga et al. (2011) case, for the corresponding progenitor parameters it is associated with the temperature of the envelope at an optical depth ∼ 1000.Blinnikov et al. (1998) arrive at the same conclusion, noting that that particular model features a very weak absorption opacity, and as a result, unreasonably high emission temperatures.
The agreement of our results with these earlier calculations pro- vides additional support to our code's validity (in particular to the applicability of the diffusion approximation), and to the conclusion of the detailed analysis of SW17, that the shock cooling emission is not sensitive to deviations of the density profile from a polytropic one.

DISCUSSION AND SUMMARY
We derived an analytic description of the deviations of shock cooling spectra from blackbody for red supergiant explosions.The approximation is given by eqs.( 56) and (57).The definitions of (and equations for) all variables appearing in eqs.( 56) and (57) are given in the Appendix.The analytic description holds from post breakout ( > 3 / ), up to H recombination ( ≈ 0.7 eV), or until the photon diffusion time through the envelope becomes shorter than the dynamical time, see eq. ( 21).The analytic expressions were calibrated against a large set of numeric MG calculations, and provide an excellent approximation to the numeric results-see § 6, figures 4 and 5.
In accordance with SW17 and Paper I, we find that the results are not sensitive to metalicity and to the ratio of core to stellar radii, in the range c / = 10 −1 − 10 −3 .In § 7.1 § 7.2 we showed the effects of deviations from and excitation LTE and of 'expansion opacity' corrections are small.Our analytic formula depends on the same four parameters as the previous one in Paper I, { , v s * , , env }.Of these parameters, the SED is most sensitive to the progenitor radius , and the ejecta velocity v s * , and insensitive to the parameters and env , where and env are the total mass of the ejecta and mass of the envelope, and is a dimensionless factor of order unity that depends on the inner envelope structure (see eq. 1).Therefore, the former two parameters are the ones most readily extracted from observations.We note also that deducing parameter values is hindered by the difficulty of determining col at early times, when the maximum observable photon frequencies (ℎ 10 eV) may be located below the thermal peak.For this regime, the frequency-dependent deviations from blackbody may prove to be an important discriminator between models of the Raleigh-Jeans part of the spectrum.
Our results are compared in § 8 to earlier numerical results including STELLA radiation transport calculation results for several non-polytropic profiles (Blinnikov et al. 1998;Tominaga et al. 2011;Kozyreva et al. 2020) obtained using MESA stellar evolution calculations.The agreement of our results with these earlier calculations provides additional support to our code's validity (in particular to the applicability of the diffusion approximation), and to the conclusion of the detailed analysis of SW17, that the shock cooling emission is not sensitive to deviations of the density profile from a polytropic one.
Based on our analysis of § 7.1, we find that the emergence of Balmer lines does not coincide with the recombination time but likely occurs earlier (approximately at ∼ 2 − 3 eV).While the recombined H fraction at this time is very low, and the effect of the lines on the SED should be small, the lines may be visible in a spectrum.The appearance of the first Balmer lines should not be associated with the recombination time 0.7 eV , when the temperature drops to 0.7 eV.
Our frequency-dependent formula agrees with the simulations' results with an RMS error of 20% (with the exception of the Wien tail, where deviations are larger but reflect a very small inaccuracy in radiation temperature).This error is somewhat larger than the uncertainty in our numeric results, ≈ 10%, as quantified by the difference between the simulations' results and the high frequency-resolution calculation (eq.61).Therefore, in order to best incorporate the uncertainty of our analytic model when comparing to observations, we suggest using a covariance matrix (as function of ( , )) of the residuals between the frequency dependent formula (eq.56 or 57) and our published simulations (see §Data Availability).
Figure1.Tests of our multigroup code.Top: an example of a radiation mediated shock (RMS) in a uniform medium, with photon energy density / ℎ , shown at several locations across the shock.This is compared with an analytic Wien profile, defined by the analytic bolometric energy density , and the local photon number density , extracted from simulation.The shock is of velocity =10%, in a pure Hydrogen medium, with initial proton number density of p = 10 15 g cm −3 .Bottom: luminosity / ℎ , at shock breakout (over the course of several shock-crossing times 0 ) with a solar composition and a free-free opacity only.The simulation is compared at breakout time to the table results of SKW III, where we extract breakout velocity, density, and breakout time based on a fit of the bolometric luminosity to SKW I (as is done in Paper I), without any additional fitting.The progenitor radius is = 10 13 cm, core and envelope masses are env = core = 10 ⊙ , and the explosion energy is = 10 52 erg cm −3 .

Figure 2 .
Figure 2. A comparison of our opacity table code results against TOPS and OP results for pure Hydrogen, showing excellent agreement.At top, a comparison of the Rosseland mean at different densities and temperatures.At bottom, example frequency dependent opacities at a density of = 10 7 g cm 3 .

Figure 4 .
Figure 4. Observed luminosity/ ℎ , shown for a range of progenitor radii and core/envelope mass combinations, comparing our analytic formulas with multigroup simulation results.Colored solid lines represent simulations.Our analytic formula from Paper I (eqs.13, 14) is in black, and our modified frequencydependent formula (eq.56) is in dashed colored lines, showing excellent agreement with the simulations.The blue / red / yellow color denotes results at times 3 / , 1.5eV , 0.7eV .In plots with with black triangles, the time 0.7 eV (in yellow) is replaced by the transparency validity time tr / , since the latter occurs prior to recombination.The formula is based on breakout parameters, ( bo , bo , ), extracted from simulation without additional fitting involved.

Figure 7 .
Figure7.Rosseland mean-free-path relative to expansion mean-free-path as a function of scattering optical depth es .Different photon frequency groups are shown in different colors, with the • and symbols indicating the location of the diffusion depth and thermal depth (estimated without expansion opacity) respectively.Results are shown at the time of recombination 0.7 eV , when the effect of expansion opacity is strongest, for two different progenitor radii, = 3 × 10 12 cm (top) and 10 14 cm (bottom), with explosion energy = 10 51 erg, and core and envelope masses core = env = 10 ⊙ .

Figure 8 .
Figure8.A comparison of the flux obtained in the numeric simulations, that do not include the effects of expansion opacity, with the flux obtained using the analytic approximations of eqs.(42) and (45) with col, (eq.23) calculated with the density and temperature profiles obtained in the simulation, but with a high resolution, Δ / ∼ 10 −5 opacity table and taking into account the Doppler shifts of the lines (The gray approximation, eqs.(13)-(15) shown in black).The high-resolution spectrum has been averaged over frequency bins, without affecting the SED.Top (bottom): progenitor radius = 3 × 10 12 cm ( = 10 14 cm), explosion energy = 10 51 erg, and core and envelope masses env = core = 10 ⊙ .

Figure 9 .
Figure 9. the photon energy distribution compared to the Planck distribution at the local plasma temperature.Photons match the local temperature deep in the ejecta and for frequencies ℎ > 3 eV.Elsewhere, they are lower than the local Planck distribution because they are dominated by photons that have arrived from deeper in the ejecta and have been suppressed by scattering.The example shown is of an explosion with progenitor radius = 10 13 cm, core and envelope masses of env = core = 10 ⊙ , and explosion energy = 10 51 erg, shown at recombination time 0.7 eV .

Figure 10 .
Figure 10.Example relaxation rates in the plasma as function of scattering optical depth es shown at recombination time, 0.7 eV , when relaxation rates are slowest.For clarity, electron collisional (photon) processes are shown in solid (dashed) lines.In the top panel, time scales for ionization and for velocity distribution thermalization by Coulomb interactions are shown.At bottom, electron collisional excitation and photoexcitation time scales are shown.Thin (thick) dashed lines denote maximum (average) photon excitation timescales based on an average of the atomic lines (see text).The characteristic time scales of all processes, except electron collisional ionization, are much shorter than the dynamical time.

Figure 12 .
Figure 12.Progenitor density profiles used in earlier shock cooling calculations, compared with the polytropic approximations used in our simulations.Tominaga et al. (2011) at top and Kozyreva et al. (2020) at bottom.The profiles we are approximating are shown in blue, and our approximation profile is shown in dashed lines.

Figure 13 .
Figure 13.A comparison of the specific shock cooling luminosity, = / , obtained by Tominaga et al. (2011) in color (different colors represent different choices for the bound-free opacity) and our MG diffusion result in black.Dashed black lines show results for opacity including only bound-free contributions, solid lines show results including both bound-free and bound-bound contributions.See § 8 for discussion.

Figure 14 .
Figure 14.A comparison of the specific luminosity obtained in earlier calculations, Kozyreva et al. (2020) (top) and Blinnikov et al. (1998) (bottom), with our results.Top: SED during shock breakout.Our results are shown for calculations with and without inelastic Compton scattering.Bottom: Shock cooling emission at different times, marked on the curves in units of days.To reduce clutter, we show our results at short wavelengths ( < 3000Å) only at 0.281 days (blue) and 3.08 days (green).For longer wavelengths, we also include intermediate times (other colors).Dashed (color) lines show results for opacity including only bound-free contributions, solid lines show results including both bound-free and bound-bound contributions.See § 8 for discussion.

Table 1 .
Observed luminosity, same as fig.4,showing example checks of the sensitivity to opacity changes.At top, comparison using our table and TOPS, showing ∼ 10% difference in col at late times, and otherwise minimal effect.At bottom, the effect of changing metallicity, exhibiting negligible effect.Both panels are showing simulations with the parameters = 10 13 cm, = 10 51 erg, and core and envelope masses of core = env = 10 ⊙ .List of lines removed from our finalized simulations (see § 7.1).

Table 2 .
Previous WorkProgenitor Radius [cm] Core Mass [ ⊙ ] Envelope Mass [ ⊙ ] Explosion Energy [erg] Parameters of the polytropic models that we use to approximate existing STELLA simulations in the literature.