Type IIP Supernova IV. Shock Breakout from Progenitor Stars Modeled with Convective Overshoot and Mass Loss

Red supergiant stars lose a lot of mass in slow winds that forms a circumstellar medium (CSM) around the star. When the star retains a substantial hydrogen envelope at the time of explosion, it displays characteristic light curves and spectra of a Type II plateau supernova (SN), e.g., the nearby SN 2013ej. When the shock wave launched deep inside the star exits the surface, it probes the CSM and scripts the history of mass loss from the star. We simulate with the STELLA code the SN radiative display resulting from shock breakout (SBO) for a set of progenitor stars. We evolved these stars with the MESA code from their main-sequence to core-collapse phase using diverse evolutionary inputs. We explore the SN display for different internal convective overshoot and compositional mixing inside the progenitor stars and two sets of mass-loss schemes, one the standard “Dutch” scheme and the other an enhanced, episodic and late mass loss. The SBO from the star shows closely time-separated double-peaked bolometric light curves for the Dutch case, as well as high-velocity ejecta with minuscule mass accelerated during SBO. The earlier of the peaks, which we call the precursor peaks, are compared with analytical expressions for SBO of stars. We also contrast the breakout flash from an optically thick CSM with that of the rarefied medium established by Dutch wind. We describe how the multigroup photon spectra of the breakout flashes differ between these cases.


Introduction
Radiative emission of a supernova (SN) 5 exploding out of a star with a large hydrogen-rich envelope has several distinct phases. While the most notable stage is the eponymous "plateau" phase for the Type II plateau (IIP for short) SNe, there is also the initial shock breakout (SBO) phase associated with a burst of bright radiation. The latter phase lasts only a fraction of a day. It takes place when a radiation-pressuredominated hydrodynamic shock wave launched deep inside the star near the iron-core boundary penetrates the outer surface of the star. This leads to an extreme brightening for a short duration (see, e.g., Klein & Chevalier 1978;Matzner & McKee 1999;Nakar & Sari 2010;Blinnikov & Tolstov 2011;Chevalier & Irwin 2011, 2021. At early times the temperature of the shells that radiate the luminosity pulse is so hot that they give rise to ionizing, X-ray, and ultraviolet radiation. A major unsolved problem of modern stellar astrophysics concerns how the diverse types of SN explosions are determined by the nature of their progenitor stars. While some broad understanding has been achieved in the past few decades on certain subclasses of SNe resulting from stripped (hydrogen, helium, or both) envelope progenitors, there is a whole landscape of properties of hydrogen-rich Type IIP SNe that explode from supergiant stars, whose properties ultimately depend on how the underlying star has evolved. Dessart & Hillier (2019) have alluded to the difficulties of inferring the masses of progenitors from the light curves of Type IIP SN even when supplemented by spectral line information that traces the rate of expansion of the spectral-line-forming regions. In particular, the discrepant estimates obtained for progenitor (and their ejecta) masses obtained via light-curve modeling of radiation hydrodynamic simulations and those from pre-explosion imaging and/or host stellar population studies are well known (see, e.g., Van Dyk et al. 2003;Utrobin & Chugai 2008;Smartt 2009;Smartt et al. 2009). The predominant reason for the similarities (in other words, degeneracies) of light curves and spectroscopic developments over timescales of 10-200 days is ascribable to the fact that progenitor stars of quite different masses at the mainsequence stage can end up retaining comparable hydrogen-rich envelopes at the time of core collapse, thus making light-curve behavior nonuniquely related to main-sequence masses of the progenitors.
One can investigate the inverse problem: what kind of radiative display of SNe results from a set of progenitors evolved, in different possible ways, from the central hydrogen-burning stages all the way to the end of the most advanced nuclear fuel burning and subsequent core collapse? This is the path we follow in this paper and its companion paper (hereafter Paper V), which investigate the light curves and spectral signatures both at the SBO and its immediate aftermath and during the subsequent phases of luminosity plateau and radioactive decay powered decaying tails. We utilize in these works our suite of progenitor models that we have evolved employing a variety of input physics and computational control that regulate the evolutionary sequence of these stars (see Das & Ray 2017;Wagle et al. 2019. This is undertaken with the Modules for Experiments in Stellar Astrophysics (MESA code) Paxton et al. (2011Paxton et al. ( , 2013Paxton et al. ( , 2015Paxton et al. ( , 2018. In the process of breakout, the radiation-pressure-dominated shock, 6 generated after the core rebounds, drives the outward expansion of the stellar envelope, and the optical depth of the material ahead of the shock decreases with time. Radiation begins to escape out of the shocked zone once this optical depth falls below c/v(r), where c and v(r) are the velocities of light and the fluid, respectively (Ohyama 1963), and breakout happens when the shock reaches the edge of the star. However, in many cases, there is evidence of a significant circumstellar medium (CSM) formed by the ejection of mass from the progenitor star prior to its explosion either in a steady wind or in episodic, massive ejections. In the cases where a dense CSM exists, the SBO may take place out of the CSM at larger distances and at later times (Chevalier & Irwin 2011, 2021Waxman & Katz 2017). Properties of this early radiation carry signatures of the progenitor star (e.g., its radius) and its mass-loss history a little before its explosion. Studying the characteristics of this early radiative emission for a given set of progenitors and their mass-loss history could therefore be of interest for the overall problem of connections between progenitors and SNe.
The final masses, radii, and structure of the pre-SN star and the energy delivered, post-bounce, at the edge of the protoneutron star affect the post-explosion dynamics and development of light curves and spectra significantly (see, e.g., . In this paper we explore how convective overshoot in the pre-SN star affects model light curves of an SN especially during the SBO phases, through scales controlled by the energy of the explosion and the extent of overshoot mixing. Higher overshoot factors generally led to more extended stars (though there are exceptions noted in , and so the underlying effect of overshoot may be reflected in both luminosity and spectral signatures at the SBO phase. We explore these results in the framework of MESA (version r-10398; Paxton et al. 2018) to initiate the core collapse and explosion and the subsequent propagation of the shock wave through nearly the entire star and then handed over to the STatic Eddington-factor Low-velocity Limit Approximation (STELLA) code (Blinnikov et al. 1998(Blinnikov et al. , 2006Blinnikov & Sorokina 2004;Baklanov et al. 2005) to complete the SBO through the star and compute radiation hydrodynamic variables in the subsequent SN stages.
We also explore the sequence of events connected with SBO and late light curves in the context of varied possibilities of mass-loss scenarios, e.g., within the "Dutch" mass-loss scenarios (see Paxton et al. 2015; Appendix A of Das & Ray 2017; Section 2.6 of Wagle & Ray 2020), as well as substantially enhanced late-stage episodic mass-loss scenarios. In particular, we find examples of high-luminosity breakout episodes that breach the classification of "superluminous SNe" (SLSNe) even if for a very brief while. In addition, we find quite generic examples of closely spaced double peaks in bolometric light curves in the "Dutch" mass-loss cases, which are usually seen in completely different types of SNe (Type IIb SNe, where the peaks are separated usually by large time differences).
Our main goal here is to explore the effects of convective overshoot and intense mass loss from the progenitor star on the variations on light-curve properties. Apart from the mass and extent of the hydrogen-rich envelope, there are other factors that affect the structure of the progenitor star, which can in turn influence the luminous display of the SN, such as the overall metallicity of the gas from which the pre-SN star formed and how convective transport inside the star is modeled. As in our previous publications (Das & Ray 2017;Wagle et al. 2019, we take the metallicity environment of the host galaxy NGC 628 of SN 2013ej as a typical situation (with Z = 0.006), which is quite similar to that of the Large Magellanic Cloud (LMC). 7 In these works on the pre-SN evolution and structure of the star, we used the publicly available code Modules for Experiments in Stellar Astrophysics (MESA) to investigate the properties of pre-SN stars and assess the dependence of stellar properties at various evolutionary stages on the convective overshoot and resultant chemical mixing inside the star. As mentioned there, the overshoot of the convective layers refers to the extra mixing of chemical elements beyond the convective boundaries defined in the traditional mixing length theory (MLT) scenario, without changing significantly the thermal structure (e.g., the entropy structure) of the layers beyond the convective boundary, but only leading to chemical mixing. In all our calculations we consider only nonrotating, spherically symmetric stars.
An interesting feature of a few cases where episodic but lowmass shells exist in the optically thin medium outside the exploding star is the existence of short-lived, high-velocity shells, which are radiatively accelerated. These are not prominent in cases where the CSM is dense and optically thick. We also find nonthermal photon spectra with bumps at higher frequencies in the examples of explosions in optically thin, low-density CSM, which are again absent from the emergent spectra of SNe in dense, optically thick CSM.
In Section 2 we recapitulate a few key attributes of the progenitors (and their hydrostatic evolution) that are important for the later explosive developments and summarize the suite of pre-SN models used in this paper. In Section 3 we outline the physics of SBO of a star with a rare and optically thin CSM, as well as when the star explodes within a dense, optically thick medium. In Section 4 we describe the methods followed to simulate the explosions and the details of how we handle the MESA and STELLA codes to generate the light curves and other radiation hydrodynamic properties. In Section 5 we present results of this study. Our results are discussed in Section 6, where we also consider different approaches adopted in the literature between gray and multigroup radiation hydrodynamics. We present our conclusions in Section 7.

The Progenitor Stars and Their Circumstellar Mediums
The final envelope mass left on the star is dependent on how much mass loss the star underwent in the stages prior to its 6 A shock front propagates with supersonic speed with respect to the material ahead of it and is subsonic with respect to the material behind it (upstream). MESA defines the shock position as the outermost location where the fluid Mach number, measured in the rest frame of the star, exceeds unity. 7 We note that for modeling Type IIP SN light curves, Goldberg et al. (2019) use pre-SN stellar models with mainly solar metallicity (Z = 0.02), except for two cases, SN 2005cs and SN 2009N, with lower-metallicity environments. Our suite of pre-SN models have a systematically subsolar (LMC-like) metallicity environment. collapse and explosion. In some cases, in the immediate past years or months, during which there could be episodic enhanced, heavier mass loss could ensue (compared to the more steady and lower rates during long-lived nuclear burning stages of the core). The lost mass carried away in a relatively slow wind (from red supergiant (RSG) stars) actually provides the CSM that influences the light-curve development of the SN when the star explodes. Unfortunately, the physical processes that control mass loss are sparsely understood, and only empirical mass-loss formulae are employed to describe the evolution of the star and the remaining stellar envelope (for a discussion of typical mass-loss rates used in our calculations and those by other groups see Morozova et al. 2016;Das & Ray 2017;. The mass loss from the star gives rise to the CSM properties that are critical inputs to light-curve modeling. Type IIn SNe (e.g., SN 2010jl) may be explosions of stars that underwent strong mass loss ( > 3 M e ) before the event and where the optical luminosity is powered by circumstellar interaction (Chevalier 2014;Fransson et al. 2014;Ofek et al. 2014;Chandra et al. 2015). SN SBO interaction with dense CSM may partly power even the Type IIP/L SNe (e.g., SN 2013ej). The dense CSM may have been due to large mass loss just (∼1 yr) before core collapse (e.g., Morozova et al. 2016;Nagy & Vinkó 2016). Note, however, that ChandraX-ray data on SN 2013ej indicate that the mass loss is well constrained on a timescale 40-400 yr before core collapse (Chakraborti et al. 2016) and has a substantially lower mass-loss rate than implied by Morozova et al. (2016) and Nagy & Vinkó (2016). The star had a mass of approximately 13 M e (Das & Ray 2017). This SN took place in the galaxy M74 (NGC 628) only about 9.6 Mpc away. Das & Ray (2017) constructed the CSM around its progenitor from the information of that episodic mass ejection during the late-stage evolution of the star in a selfconsistent way. The SN belonged to a subclass identified as a moderately interacting SN intermediate in class between Type IIP/L and Type IIn SNe that have dense circumstellar medium immediately outside the progenitors resulting from a phase of late-stage heavy mass loss. The presence of such a dense medium based on the early light curves of these special Type IIP/L SNe was invoked by several authors (Valenti et al. 2015;Morozova et al. 2017;Yaron et al. 2017).
We summarize the properties of the progenitor stars used for this study, at the core-collapse stage, with dense CSM fed by enhanced mass loss and CSM generated by Dutch wind mass loss in Tables 1 and 2, respectively.

Shock Breakout from the Star and Its Circumstellar Medium
To contextualize our results, we briefly review here the literature on the physics of SBO from the progenitor star with a rarefied (optically thin) CSM or from a dense and optically thick CSM. While in the former case the SBO can be considered in a plane-parallel geometry (at least for short times immediately after SBO), in the latter case, where an extended dense wind is present, SBO is analyzed in spherical geometry (Balberg & Loeb 2011;Ginzburg & Balberg 2012). We first review the case of SBO of a bare star or that with a low-density CSM. Then, in Section 3.3 we discuss the case of SBO of a dense, optically thick CSM created by previous episode(s) of heavy mass loss from the progenitor star sometime before the collapse and explosion.

Self-similar Solutions of Shock Propagation and Acceleration
Chevalier (1992), Chevalier & Soker (1989), and Imshennik & Nadezhin (1989) had used the self-similar solutions of Gandel 'Man & Frank-Kamenetskii (1956) and Sakurai (1960) to describe the shock acceleration in a power-law density profile as a function of the radial coordinate from the edge of the star to describe the development of the asymptotic power-law density distribution after the SN explosion 8 (see Equation (2.14) of In this way, these authors describe the outer atmosphere of the exploding SN 1987A resulting from a blue supergiant, in the form of a power-law density profile. Type IIP SNe arise out of RSG stars (Smartt et al. 2004). Matzner & McKee (1999) describe the relation between the structure of the progenitor star and the distribution of the ejecta set in motion by the effect of the shock wave that initiates the SN explosion. 9 This relation is established by the hydrodynamics of the explosion, which consists of several separate stages. First, the shock wave propagates through the stellar mantle and envelope, compresses this material, and sets it in motion, a stage called the blast wave, which propagates initially in an optically thick medium (Chevalier 1976). The blast wave, a radiation-pressure-dominated strong shock, in the medium with a power-law density profile (progenitor structure approximated by ρ = Ar − w ), is described by adiabatic, self-similar solutions (Taylor 1950;Sedov 1959). The character of the solution varies with the progenitor density index w itself (typically w ∼ 1.5-3 above the core). The progenitor of the SN prior to the collapse has different composition zones in which the density and radial scales vary widely. The blast wave going through each layer decelerates inside each layer. However, in between compositional layers (like the gaps in onion-skin structures) there are steep density gradients between zones. The shock accelerates while passing through these gradients-this is illustrated in Figure 2 of Matzner & McKee (1999). At the end of the blast wave phase, when the shock approaches the surface of the star, it rapidly accelerates. Chevalier (1976) considered the structure of the blast wave in an increasingly steep (progenitor) density distribution and likened this to a progression of Sedov blast waves in different power-law mediums. Equation (19) (or, equivalently, Equation (17)) of Matzner & McKee (1999) provides an approximation to the shock velocity v s in terms of the initial, local Lagrangian variable, the mass m of the enclosed ejecta, initial radius r 0 , and density ρ 0 , as well as the energy E in of explosion. The shock velocity was obtained by an interpolation between Sedov and Sakurai selfsimilar solutions and with the assumption that the opacity κ is not space and time dependent (see also Rabinak & Waxman 2011, hereafter RW11, especially their Equation (3) 10 ).

Shock Breakout Pulse and Deviations from Thermal Equilibrium
The shock wave propagates as an adiabatic shock, i.e., without energy losses, as long as the optical depth of the material ahead of the shock exceeds a certain critical value. The shock wave is a radiation hydrodynamic structure in which photons can diffuse upstream and get amplified in energy by compression. For this to continue to happen repeatedly, the time for light to diffuse across the shock front has to be at most comparable to the time for the fluid to flow through it, in other words, for as long as the optical depth for photons to escape the remaining stellar envelope τ(r) is larger than c/v sh , where c is the velocity of light and v sh is the velocity of the shock at radius r. When the overlying optical depth of the shock τ s becomes comparable to ∼ c/v sh , the shock is said to break out and leave behind an expanding sphere of gas and radiation (Ohyama 1963;Chevalier 1976;Weaver 1976;Matzner & McKee 1999;Nakar & Sari 2010). The shock usually has high velocities v sh /c 0.01 at these times going through ionized plasma at densities ρ ∼ 10 −10 to 10 −6 g cm −3 . Under these conditions, free-free and bound−free emissions generate photons, and the opacity is largely due to electron scattering (Thomson scattering). The kinetic energy of the gas and radiation mixture is converted to thermal energy via Compton scattering of photons on free electrons (Sapir et al. 2011). Nakar & Sari (2010) argue that when the luminosity shell is deep inside a highly optically thick region (  t 1 ), the The distribution of radioactive matter in the inner ejecta may also be a complicating consideration at a slightly later time. 10 Caveats mentioned by Ganot et al. (2020) about RW11, with the modifications of Sapir & Waxman (2017, hereafter SW17), state that RW11 describe the shock cooling phase of an SN and do not treat the interaction between the SN radiation and a circumstellar medium (CSM) and therefore can be applied only for SNe with no evidence for CSM interaction in their spectra.
escaping photons would have undergone many interactions with electrons as they travel across the envelope until they are able to escape after reaching a layer where the ambient τ ∼ 1. Thus, the typical energy of each of the escaping photons is controlled by photon-electron interactions along their way outward from the luminosity shell. The important processes that control whether the radiation is in thermal equilibrium with matter in the outflow are free-free and bound−free emission and absorption, as well as Compton and inverse Compton scattering. They define the thermal equilibrium temperature T BB in terms of a photon energy density for mass shells external to the luminosity shell (which has an optical depth t and luminosity L) up until the optical depth falls to unity, which is with "a" being the radiation constant. Nakar & Sari (2010) assume that the free-free process dominates the process by which an isolated shell with radiation-dominated energy density ò tries to approach thermal equilibrium. The time required to achieve thermal equilibrium in the isolated shell is approximately  n n T BB ph,ff BB ( ), where the denominator is the rate of generation per unit volume via free-free emission of photons with energy hν ∼ 3kT BB . Thus, a thermal coupling coefficient η for the expanding gas is defined by Nakar & Sari (2010), in terms of the available time t t min , d { }(where t d is the diffusion timescale) that photons are confined to the particular isolated shell. However, the same group (Shussman et al. 2016) modified this equation (Equation 9 -see the footnote 4 of Shussman et al. 2016) with an overall numerical coefficient on the right-hand side that is 4 times higher, due to the inclusion of bound−free transitions in opacity (in addition to the free-free opacity), resulting in a higher breakout shell temperature, and a (much) shorter photon diffusion timescale. The changed result is Here T obs,0 is the typical photon energy of the observed radiation, but its spectrum is not necessarily a blackbody. The above definition of η 0 , however, does not include Compton boosting of low-energy photons, which is important when matter and radiation are out of equilibrium. For systems of interest here, where T BB 100 eV, Comptonization is important and should be included in determining the radiation temperature T obs once the radiation is out of equilibrium, a la Equations (11)-(15) of Nakar & Sari (2010; see also Waxman & Katz 2017). On the other hand, if η 0 > 1, the shock at breakout is not in thermal equilibrium. Given the weak dependence of η 0 on density ρ 0 , this condition for thermal equilibrium translates to the shock velocity at breakout v 0 < 11,000 km s −1 . As we shall see in Section 5, the breakout flash is well out of thermal equilibrium according to this criterion. This is also evident from the spectral shape of the emitted radiation (see Figure 12). Nakar & Sari (2010) point out that when η > 1, the spectrum is an optically thin thermal freefree emission (with electrons at temperature T) modified by Compton processes and consists of a Wien spectrum at higher frequency and a thermal spectrum of temperature T at low frequency. The spectral shape of the emitted radiation (see Figure 12) is actually complex, showing a separate high-energy component due to photons from hotter layers becoming visible at breakout stages.
Between the planar and the spherical phases, the nature of the radiation changes. During the planar phase, the optical depth of each mass shell remains nearly constant, while diffusion time from each mass element grows linearly with time, as does the dynamical time R/v sh also, and thus radiation escapes from the same mass shell from which the radiation initially escaped at the time of breakout. On the other hand, during the spherical phase, the optical depth of each mass shell keeps decreasing with time because of the spherical dilution factor, and therefore deeper and deeper layers start dominating the emitted radiation. We adopt the definition t s = R * /v 0 .
The typical timescale associated with the SBO from an RSG star is the same as the diffusion time at breakout t 0 given above. At the early part of the rise of luminosity, some photons diffuse ahead of the shock, even when the shock is somewhat far away from the edge at x = − vt. As the fraction of photons that diffuse a length x ahead of the shock is proportional toe x 2 , the early rise is argued to contain a terme t 2 . The radiative energy itself around breakout phase scales as e − x . The shape of the light curve is thus represented by (Sapir et al. 2011;Shussman et al. 2016 where L 0 ≡ E 0 /t 0 and E 0 is the internal energy of the breakout shell; see Equations (34) and (40) of Shussman et al. (2016) for expressions of L(t) and thus L 0 . They also obtain a best fit to their numerical results for p = 0.35 and q = −0.15. Our approximations to the numerical results of luminosity pulse obtained by STELLA models by similar functions are discussed in Section 5.2. Calzavara & Matzner (2004) point out that there is a correlation between the photon temperature T se and diffusion time t se (equivalent to t 0 in the notation of Nakar & Sari 2010), which is also termed the time for shock emergence. The temperature that Calzavara & Matzner (2004) refer to is equivalent to T BB (at the breakout shell) and T rad in outputs of the STELLA code and in the notation of Nakar & Sari (2010) and Shussman et al. (2016). This is given by Here a = 4σ/c is the radiation constant. For a given progenitor (i.e., with a specified constant "n"), all the terms on the righthand side are essentially constants (with an appropriate κ assumed by Calzavara & Matzner 2004 to be constant for the given metallicity of the progenitor environment).

Type IIP SN Shock Breakout in a Dense CSM
In the previous section we outlined the process of SBO of a bare star and the consequent emission of light pulse and its thermal character (or lack thereof). For an RSG star, given its large size, it takes a good part of a day for the shock to travel through the star to its edge. The duration of the SBO and the light pulse immediately after is much shorter-of the order of 10 3 s (Klein & Chevalier 1978). However, during the long evolutionary process of the progenitor, the star may have lost mass in a wind that, because of its relatively slow speed, still makes up a surrounding medium (the CSM) that is denser than the normal interstellar medium. If the CSM density is low, then the CSM, through which the shock eventually travels after it exits the star, is optically thin. In such cases, many of the characteristics of the SBO and its accompanying light pulse (outlined in the previous subsections) remain qualitatively similar, although there may be a few unique properties, like the ejection of a radiatively accelerated highvelocity shell of minuscule mass for the Dutch mass-loss case (see Section 5.4). In other cases involving Type IIn and some specific subclasses of Type IIP SNe, the progenitor star may have undergone late-stage, intense mass loss, just before its core collapse, leading to a dense CSM. When the shock is sufficiently far into the CSM as the overlying optical depth of the CSM becomes comparable to c/v sh , radiation begins to leak out of the CSM. The dynamics and bolometric light curves of SBO from a dense CSM (which are qualitatively different from those of a bare star or of a star with optically thin CSM) were studied using analytical models (Balberg & Loeb 2011;Chevalier & Irwin 2011, 2021Chatzopoulos et al. 2012Chatzopoulos et al. , 2013 and numerically by Moriya et al. (2011) and Moriya (2013), among others. For the high-density CSM, the gas, initially, is most likely neutral and dusty. The dust opacity may be the dominant opacity when dust is present. But the radiation-dominated shock wave has a precursor in the CSM region that is expected to heat and evaporate this dust soon and reduce the opacity. Furthermore, the photosphere is bounded by the region where there is a sharp gradient in opacity, as the gas is also ionized. This ionization front also moves out as the photosphere moves outward through the dense CSM. Once the CSM is ionized, the motion of the photosphere slows down and the temperature rises. After this, the luminosity and temperature rise sharply until their maxima. Chevalier & Irwin (2011) distinguish two cases for the radiation-dominated SBO from an opaque CSM based on the relative magnitudes of the spatial scales R w (the outer boundary of wind-fed CSM) versus R d , the characteristic diffusion radius (or the radial distance the shock has to travel after exiting the star from which leakage of radiation starts). This R d is obtained for a specific CSM density stratification, namely, for a steady wind with a constant wind speed (v w ), . For the dense massloss case, the time taken by the shock to move to the breakout region, or the onset of light flash, is The rise time for the light curve has been shown (Chevalier & Irwin 2011) The region around the SBO is not expected to be extended in radius, and the radiation pressure of the escaping radiation accelerates gas until the radius R w . A dense shell may form where the radiation escapes, but the dense gas does not extend much beyond the SBO radius, which is close to the external boundary of the wind. The pulse shape schematically drawn in Figure 1 of Chevalier & Irwin (2011) shows the pulse width to be roughly twice the rise time. The (corrected) rise time and the luminosity peak at breakout are given as functions of ejecta mass, energy of explosion, and other parameters as 0.34 cm g 10 days Here D * is related to the normalization D of the CSM density profile by D * = 0.05D/c. Note that the product L Peak × t rise 1/2 is independent of D * and has the strongest variation with the energy of explosion:

Methods and Key Aspects of the Simulations
The starting point of our numerical simulations with MESA and STELLA are the progenitor star structure and compositional variables determined earlier until the onset of core collapse (Wagle et al. 2019. Our chosen progenitor stars are nonrotating stars with metallicity Z = 0.006 starting from zero-age main-sequence (ZAMS) masses in the lower range of massive stars (12-14 M e ), consistent with the expectations for the progenitors of the prototypical Type IIP SN SN 2013ej in this context (see . Briefly, these progenitor models with no rotation were evolved with MESA using the soft-wired 79 isotopes network (Farmer et al. 2016). Convection was treated in the MLT approximation with α MLT = 2, we used the Ledoux criterion for determining convective boundaries, and the semiconvection efficiency parameter was α SC = 0.1. The final masses and structures at the corecollapse stage are diverse, depending on the input physics and parameters fed to MESA. These are primarily the different massloss scenarios, especially the late-phase intense mass-loss episodes and how convective transport and compositional mixing take place inside the star. For the mass-loss scenarios in the cases where enhanced mass loss is skipped, we choose the wind scaling factor η Dutch = 0.5 for the "standard" Dutch mass-loss scenario described in Section 2.6 of .
In addition to the models with standard mass-loss rate discussed above, we considered a grid of models in the same range of initial masses and overshoot parameters as above, but with an enhanced mass loss in the late stages of evolution. Das & Ray (2017) demonstrated through simulations of explosion of a 13 M e star that the light curves of SN 2013ej are better fitted in the case where they considered enhanced mass loss at the end stage of the evolution. The enhanced mass loss could be triggered by wave heating during core Ne or O burning phases. Fuller (2017) shows through simulations of a 15 M e , nonrotating, solar-metallicity star that a small amount of mass (1.0 M e ) may be ejected with low velocity (50 km s −1 ) around a year before core collapse because of wave heating during the core Ne burning phase. In our models, we implement such an enhanced mass loss during the core O burning stage, by setting a condition on central mass fractions of oxygen and sulfur through run_star_extras. Enhanced mass loss of 0.7 M e yr −1 is applied after center O burning begins, provided that the mass fraction of 32 S at the center is >0.03 and the mass fraction of 16 O at the center is >10 −3 . This condition triggers an enhanced mass loss of 0.7 M e yr −1 during the last 2-4 yr before collapse, for our grid of models. For these enhanced mass-loss models, we used the 22 isotopes network (see Section 2.4 of Wagle & Ray 2020, hereafter Paper II).
We specify the various overshoot factors f used in our pre-SN models constructed in our previous publications (other convection-related parameters f 0 and mixing length parameter α MLT are kept fixed as in our previously published stellar models). The properties of the progenitor stars are given in Tables 1 and 2. The numerical modeling of collapse and explosion is carried out in two steps. First, the infall, explosion initiation, and excising of the core, followed by the shock propagation throughout the star to a zone close to the edge of the star, are handled by MESA (as described in Section 6.1 of Paxton et al. 2018). We study the effects of varying amounts of energy injected at the shock generation stage on the explosion development, including the SBO, by changing the relevant parameters in the inlist file: inlist_edep (see Section 4.1). In the second stage, once the shock has reached a designated mass shell close to the edge of the star, the computational control is handed over to the STELLA code, which develops the hydrodynamic evolution of the star through the SBO phase and the subsequent evolution as an SN in its plateau phase until well into the radioactive decay phase of its light curve.
We use as templates the inlist files provided in the test suits package in MESA such as the example_ccsn_IIp (this appears in the r-10398 release package subdirectory star > test_suite). Unless stated otherwise, we use the default versions of these inlists. We keep the amount of Ni 56 injected in the mantle region fixed as 0.028 M e .

Initiation and Propagation of the Shock Wave
Since MESA cannot model the transition from the core infall stage to the explosion itself, the code provides an approximate way of generating an explosion via a shock wave by injecting energy in a thin layer (0.01 M e ) and in a very short time (5 ms), roughly near the boundary of the infalling core (see Section 6.1 of Paxton et al. 2018), which has a compositional transition from Fe to Si/S. This is the so-called "thermal bomb" scenario (Aufderheide et al. 1991;Bersten et al. 2011;Morozova et al. 2015). Our choice of the duration of energy deposition is close to that of Morozova et al. (2016). The thickness of the deposition layer chosen here is similar to that in Morozova et al. (2015). The strength of the explosion is varied by setting the parameter inject_until_reach_model_with_to-tal_energy in inlist_edep to the desired value. We take a range of explosion energies (EDEP) typically between about 0.6 and 1.4 FOE.
After the completion of the energy deposition step, MESA propagates the shock wave through the star across several compositional barriers and associated steep density gradients. These are executed with inlist_shock_part1 through inlist_shock_part5 with its default parameters in the mesa/star/test_suite/example_ccsn_IIP directory of the MESA release v-10398. During these phases, the shock traverses an enormous range of density and temperature inside the star. When the shock encounters steep density gradients, such as in the He core/hydrogen envelope boundary, it splits into a forward shock/reverse shock structure. Not only does the reverse shock manage to decelerate the expansion of the core, but its concomitant increase of entropy (and temperature) also heats up the inner material, compensating the cooling effect from the overall expansion due to shock passage (Paxton et al. 2015).

MESA to STELLA Handoff
When the shock has traversed nearly the entire H envelope and is near breakout, we allow MESA to transfer the control of numerical evolution of the shock to another code, STELLA, as described in Section 6.3 of Paxton et al. (2018). STELLA handles the SBO more accurately, being a multigroup radiation transfer code, than the gray MESA code at or after the breakout, especially if any matter is placed above the photosphere (as is the case of a CSM) or when substantial radiation starts free streaming from just below the photosphere even prior to SBO. We chose the handover point to STELLA as when the outgoing shock reached a mass coordinate measured inward from the surface of the star, δm = 0.05 M e , by setting the parameter x_ctrl(16) in inlist_shock_part_5. STELLA does not assume homologous expansion of ejecta, and so early handovers can be properly executed.
The handover is facilitated by creating two input files for STELLA, which are output files of MESA, namely, mesa. hyd and mesa.abn. The former file contains information about the structure that MESA wishes to pass on to STELLA, while the latter file contains information about the composition. These files are created in the mesa/star directory by setting controls in &star_job.

Radial Profile of the Circumstellar Medium
For each model, a customized density profile for the circumstellar mass (CSM) was generated using the data from the evolution run in MESA. The mdot parameter gives the total mass of material ejected from the star in one step, while the dt parameter gives the duration of the step. Using this and assuming a constant wind speed of 10 km s −1 as the rate at which the ejected material moves outward (starting from the stellar surface), a time evolution of the CSM density profile is performed. This is necessary in order to account for the changes in stellar radius, which may at times occur faster than the wind speed. If the change in stellar radius is negative (surface of the star moves inward), the model does nothing, but if it is positive (stellar surface moving outward), the model assumes that the entire CSM present at the time is pushed outward along with the star. By this method, the innermost layer of the CSM cannot, at any time, be at a radius smaller than that of the stellar surface.
The density profile is constructed in each case taking the mass loss from the last 10 yr of evolution, which gives a density profile up to a distance of ≈5000 R e from the stellar center in the case of the Dutch mass-loss scheme. This is a sufficient extent for our modeling purposes, consistent with radii for extended CSM breakouts expressed in Waxman & Katz (2017) and in the lower range of dense wind outer radius expressed in Moriya (2013).
For the models with the enhanced mass-loss scheme, the section of CSM formed by the periods of enhanced mass loss has density values many orders of magnitude higher than the rest of the CSM. For these models, we introduce a cutoff, restricting the CSM profile to the part with density greater than 10 −14 g cm −3 in order to avoid numerical difficulties in the STELLA runs.
The completed profile has of the order of 10,000 data points. This is then uniformly sampled over radius (by linear interpolation) to match the number of lines allocated for the CSM in the STELLA input files. This is sufficiently accurate since we are sampling of the order of 1000 points from a much higher number. Uniform sampling is necessary, as STELLA has numerical difficulties and fails to converge even after weeks with nonuniform sampling. The values for density, radius, mass coordinate, and mass are then written into the mesa.abn and mesa.hyd files for input to STELLA. The values of all other parameters in these files for regions in the CSM are initialized with identical values to the last zone describing the stellar interior (zone 270 in our models), and the relaxation to their actual values occurs in the initial time steps of STELLA calculations.

STELLA: Multigroup Radiation Hydrodynamics for SN Evolution
A customized interface between MESA and a public version of STELLA (inputs, files, controls, etc.) has been provided in version r-10398 of MESA. STELLA code description has been given in Paxton et al. (2018) and in the README file in the zip-file of this release. STELLA is a 1D Lagrangian, implicitly differenced hydro code with a comoving grid, which models time-dependent, nonequilibrium, and multi-frequency radiative transfer inside the SN ejecta as it moves out under the action of the shock, before the expansion is homologous. The STELLA code is thus well tuned to handle the SBO phase of the explosion.
The basic framework for the comoving frame radiative transfer and hydrodynamics is outlined in Castor (1972Castor ( , 2004 and Mihalas & Mihalas (1984). To connect to our simulations, we note, very briefly, that the STELLA code solves nonequilibrium radiative transfer equations (Blinnikov & Bartunov 1993) correct to the order u/c using angular moments of the intensity averaged over each of its 40 frequency groups. Important questions of input physics, namely, the equation of state of stellar matter including its ionization and recombination status, its opacity and its appropriate frequency averages, radioactive energy input, are discussed in this publication. All spatial and frequency derivatives in the differential equations of comoving radiation hydrodynamics are replaced by finite differences for numerical solutions. For each Lagrangian zone, the evolution of radial coordinate (r), velocity (u), and temperature (T) and the angular moments of the invariant photon distribution function 11 f (r, μ, ν) are tracked. The solution scheme is an implicit one for both radiation and hydrodynamics variables. It has an automatic choice for the time step and allows the code to overcome the Courant restriction.
The Eddington factor f E = f E (r, ν), which is required to close the equations of radiation hydrodynamics, is estimated by solving the time-independent radiative transfer equations by introducing "impact parameters" p, in the (p, z) mesh system (instead of (r, μ) mesh), and applying the Feautrier method (Feautrier 1964; see Section 2.5.2.2 of Moriya 2013 andBartunov 1993 for the appropriate outer boundary condition). The opacity is computed by taking into account over 153,000 spectral lines (Kurucz & Bell 1995;Verner et al. 1996). The line opacities are calculated taking into account the effects of high-velocity gradients. Opacity includes the usual components like free-free absorption, photoionization, and electron scattering. For the matter component, the ionization and level populations of the atoms (of 16 elements) in the plasma are obtained via the Saha equation. Radioactive energy input from Ni and Co (into positrons and γ-rays) is tracked in a one-group transport approximation.

Output Files of MESA and STELLA and Time Parameter
We predominantly use the output files mesa.lbol, mesa. tt, and mesa.swd from STELLA for SBO analyses. Available documentation on STELLA mentions that several variations of the time parameter are in use across the several output files. The mesa.swd file, which we use to generate radial profiles of variables at several instants before and after breakout, reports time in the retarded form to account for the angular spread in the luminosity shell with respect to the line of sight to the observer. The time variable is t comoving − R out /c in mesa.swd. We derive luminosity and temperature profiles from mesa.lbol and mesa.tt, respectively. Although we do not find in available documentation any explicit details on the time variable employed in these files, we observe that the instants at which data are reported in mesa.lbol and mesa. tt do agree, and hence we conclude that these files employ the same time variable. We also conclude that the times reported in these files are retarded times by noting that the luminosity of the outermost shell of the CSM at various times in mesa.swd is also the bolometric luminosity at the various times reported in mesa.lbol, up to differences due to different temporal resolution across these files. We note that, according to Kozyreva et al. (2020), STELLA does not ordinarily account for the smearing effect of the radiation pulse.
The mesa.swd.ph file gives details about the photosphere, but the time variable here is relative to maximum bolometric luminosity. In order to get retarded time from this variable, we add to this the time for maximum L bol evaluated from mesa. lbol. In our plots and discussion, we denote the retarded time in STELLA as Time in STELLA or t STELLA , which is used to compare data across several output files.
The MESA inlists data used for our simulations presented here and some of the STELLA output files are made available publicly online at https://doi.org/10.5281/zenodo.6665072.

Results
In this section we discuss first the general characteristics of the light curves resulting from numerical simulations of SBO in SNe with STELLA. This is followed by a description of the density ramp at the edge of the pre-SN star just before SBO (Section 5.1). In Section 5.2 we discuss light-curve features for both explosions inside a rare, optically thin CSM, as well as those inside an optically thick CSM. In Section 5.3 we compare characteristics of the numerical simulation of the explosion SBO inside rare, optically thin CSM with those from analytic models. This is followed by the analysis of high-velocity ejection of a low-mass shell in the CSM in the same case of an explosion inside an optically thin CSM in Section 5.4. In Section 5.5 we discuss the early-time evolution of the profiles of radiation hydrodynamic variables inside the star and the CSM as the shock propagates out. We divide this subsection into two parts: one for explosion inside an optically thick CSM (Section 5.5.1) and the second for explosion inside an optically thin CSM profile established by a prior Dutch wind from the progenitor (Section 5.5.2). In Section 5.6 we compare our numerical results with the analytic formulation of SBO inside an optically thick CSM. In Section 5.7 we describe the photon spectral energy distribution of SBO for both optically thick and optically thin CSM.
Our explosion models are relevant mainly for the planar phase of the SBO until the onset of the "spherical" phase for the Dutch cases. For the enhanced mass-loss cases, the breakout extends to the "spherical" phases. However, in both sets of models, the post-SBO evolutions are shown well before their plateau phase. In this study, we hold the mass of 56 Ni to be a constant value of 0.028 M e . This radioactive element is distributed between the excised core boundary and the helium core boundary. In early light curves the influence of radioactive decay is somewhat marginal, but it becomes critical at the end of the plateau stage. We shall investigate the effects of variation of 56 Ni mass on the light curves and spectral line development in a subsequent paper.
The bolometric light curves for the Dutch and enhanced  m cases are substantially different: not only is the peak of the luminosity delayed in the case of the enhanced  m case with respect to the similar Dutch case (with the same energy of explosion), but the former are more luminous at the peak of L bol . Moreover, while the luminosities in the Dutch and the enhanced  m cases essentially overlap with each other by the time the plateau and the subsequent radioactivity-powered phases ensue, the luminosities are substantially different for the initial 2-week period after explosion. Thus, the mass loss from the progenitor and the shock propagation inside a rarefied versus dense CSM lead to discernibly different bolometric light curves in the initial 2-week period.

Density Ramp at the Edge of the Presupernova Star at Core Collapse and Just before Shock Breakout-Dutch Cases
We display profiles of density as a function of radius from the surface for some of the Dutch mass-loss models with EDEP = 1.0 FOE, analogous to Figure Figure 1(a), which displays these profiles at the core-collapse stage, we see an unexpected yet small bump in density at log ρ ≈ [−8.30, −8.45], which is also visible in Figure 2 of Morozova et al. (2016) at about the same density. We observe that the density profiles even at the time of the core-collapse stage are more sensitive to the overshoot parameter than progenitor mass. For the 14 M e models, we note that progenitors with larger f display lower densities at core collapse. We also note from Table 2 that for the overshoot parameters explored in this figure R cc is a proxy for f. This implies that progenitors with larger core overshoot tend to have a larger radius and hence smaller surface density at core collapse.
These density profiles are substantially modified by handover time to STELLA, as shown in Figure 1(b). However, the smaller bumps are observed at almost the same density, as they are yet unaffected by shock propagation from the core to the surface. For the 14 M e , f = 0.040 star, the small density bump had smoothened out, and its entire density profile appears different from other similar cases observed at the core-collapse stage. We also observed from our analysis of possible correlations between precursor-peak luminosities and their FWHMs in Figure 5(b) that the f = 0.040, 14 M e case represents some structural difference that stands out from the rest of our runs. Here we infer that such differences possibly originate in the duration between core collapse and shock propagation to the surface. Additionally, we note that the position of the luminosity shell marked by the τ = 1.2 · c/v surface coincides with shock position marked by the outermost zone where the fluid Mach number exceeds unity at handover time to STELLA.

Early Light Curves and Their Features
In Figure 2(a), we plot L bol , T eff , and T BB as a function of STELLA time (see Section 4.5) for our Dutch mass-loss model with progenitor mass 12 M e , f = 0.025, and EDEP = 1.0 FOE. We observe multiple peaks in L bol after 0 days for all of the Dutch mass-loss models. Henceforth, we refer to the first peak in L bol , T eff as the precursor peak and to the last major peak as the main peak. We occasionally refer to the small peaks that occur between the precursor and main peaks as the intermediate peak(s). Note that the precursor and main peaks are well separated compared to the light-travel time ("retardation" time), which is typically 1400 s or 0.016 days.
We consistently take the precursor peak to correspond to SBO since the trends of its parameters are in broad agreement with the analytical expectation for the luminosity peak versus its duration of the breakout flash in Calzavara & Matzner (2004). We also find that for each of the L bol peaks there is a corresponding T eff and T BB peak at about the same time. We note that T BB is roughly an order of magnitude higher than T eff by the time of the main peak, while the difference is smaller at the luminosity trough immediately after the intermediate peak.
This indicates that during SBO, and during the subsequent luminosity peak that we identify with the interaction of the shock with the ambient CSM, the radiation is out of equilibrium, i.e., η 0 ? 1.
On comparing the precursor and the main peaks, we observe that the latter is greater than the former for the L bol peaks, while this trend is reversed for the T eff peaks. Henceforth, we shall refer to the precursor peak in L bol as L 1 and the main peak in L bol as L 2 .
Next, we extend similar analysis to the closest enhanced mass-loss model, which corresponds to 12 M e , f = 0.020, and EDEP = 1.0 FOE in Figure 2(b). Here we observe only one peak in L bol , T eff , and T BB whose peak values coincide well at 3.92 days.
We now consider the variations in L 1 and L 2 among all the Dutch mass-loss models. Figure 3(a) is a plot of these luminosities as a function of the STELLA times at which these peaks occur. We denote L 1 (precursor peaks) by crosses and L 2 (main peaks) by circles and triangles to differentiate models based on their progenitor mass. We also use a color scheme to distinguish among a subset of the L 2 points for different families of EDEP, with darker shades denoting higher energy. Interestingly, the L 1 peaks, which correspond to our identification of the breakout flash, fall in a nearly linear curve, with negative slope indicating that the earlier peaks are also the more luminous ones. At large shock emergence times, the curve tends to have a smaller slope. We also depict a subset of the L 2 points to highlight the trends in their scatter with variation in EDEP and progenitor mass. We observe that the L 2 points cluster into different bands of constant EDEP, with the larger EDEP points occupying the higher end of the luminosity axes. We show later that the L 2 peaks occur owing to shock −matter interactions in the CSM. Since higher-energy explosions give rise to more energetic shock waves, shock−CSM interactions give rise to more luminous L 2 peaks for higher EDEP. However, when displayed along with closely varying EDEP points, these clusters are not clearly distinguishable from neighboring clusters of a different EDEP value, due to the vertical scatter among the points in a cluster. Each such cluster of points with a fixed EDEP value displays a positive correlation with STELLA time at which these peaks occur, unlike the L 1 points. Relabeling the L 2 points for different families of f did not display any discernible trend.
In order to explore the inverse dependence of L1 peaks on STELLA time, we display a subset of the Dutch mass-loss precursor peaks with progenitor mass 14 M e in Figure 3(b). Each point in this figure represents the EDEP and f of that particular model run with the help of a shape and a color scheme, respectively. By focusing on points of a particular shape, we infer that the effect of varying f alone is that peak luminosities are larger and emerge earlier for larger values of f. Similarly, we observe, by focusing on points of a particular color, that the effect of varying energy alone is that higher EDEP values give rise to larger precursor peaks, which are also attained earlier in STELLA time.
Increasing EDEP gives rise to more energetic shock waves. This, in addition to producing brighter L 1 peaks, increases shock velocity and decreases the time to SBO, hence giving rise to earlier peaks. We noted from Figure 2(a) that T eff peaks at about the same STELLA time as L bol . We observe from Table 2 that for the models under consideration the overshoot factor varies monotonously with R cc , and f is thus a suitable proxy for R cc . Therefore, the precursor peaks in L bol scale proportionally with f as in Figure 3(b). We pointed out in  that the luminosity of the SN at the plateau stage scales as the radius at SBO. Here we find that even at the SBO stage the peak luminosity scales with R cc .
STELLA does not define the handover time from MESA to STELLA clearly. Because of the ambiguity of the zero of time in the STELLA output with respect to the shock initiation time at or near the Fe core edge, we use an alternate temporal measure of the SBO peaks. Instead of the times reported in mesa.tt, we use the pulse FWHM of the L bol peaks. We fit a Gaussian curve to the L 1 and L 2 peaks separately, using the scipy.optimize.curve_fit function in Python. We fit for optimal values of the parameters: mean, standard deviation and peak value of the Gaussian such that the sum of the squared residuals of (Gaussian fit) − (STELLA L bol values) is minimized. We then calculate FWHM from the fit parameter for standard deviation. However, this method is subject to the following caveats. FWHM evaluated from Gaussian fits to the Dutch mass-loss models will be smaller than a read-off-thegraph method since the Gaussian falls more sharply than the STELLA light curves after the first L bol peak. The broadening of our precursor pulse, in the region after the peak, can be attributed to heating of the CSM, which partially insulates the escaping radiation. This post-peak broadening also introduces a temporary plateau-like flattening to the light curves, before the rise to the intermediate peak. Ignoring this plateau region in the Gaussian fits gives us a better correlation between L 1 and FWHM. The subsequent analysis of the L bol peaks uses FWHM instead of STELLA time.
In Figure 5(b) we display L 1 peaks of 14 M e , Dutch massloss models as a function of their FWHM for different families of EDEP and f, distinguished by a shape scheme and a color scheme, respectively. We return to the solid-line fits to Equation (6) of Calzavara & Matzner (2004) later. For each family of overshoot parameter, we observe a strong correlation between L 1 and FWHM, similar to Figure 3(a). Within any f = constant band, we observe that larger EDEP explosions give rise to smaller pulse widths and larger precursor peaks, with the exception of f = 0.040. While the effect of varying EDEP on L 1 is the same as that observed in Figure 3(b), the effect on the relevant timescales is more evident here. Not only does larger EDEP give earlier breakouts, but it also gives shorter pulse widths by giving larger ejecta velocities.
As we move across the different f = constant bands in Figure 5(b), we observe that larger f models give rise to larger precursor peaks and broader pulses. While the effect of f on L 1 peaks can be understood from the discussion of Figure 3(b), we  Table 2 that for models that differ only by f, larger f values give rise to a larger pre-SN radius (R cc ). It is reasonable to expect stars with larger R cc to have large radius at breakout. We inferred in our discussion of Figure 2 that L bol depends predominantly only on the radius, at breakout time. For models with very large overshoot parameters, like f = 0.040, not only is the breakout pulse exceptionally broad, but the correlation with EDEP is also reversed. This indicates some significant structural differences in the pre-SN evolution of progenitors with large overshoot parameters.
A similar analysis involving the main (L 2 ) peaks displayed weaker correlations for variations in EDEP and f, thereby indicating that the physics behind the precursor (L 1 ) and main (L 2 ) peaks are significantly different.
We extend a similar analysis to the enhanced mass-loss models. We display peak luminosities as a function of their FWHM for different families of progenitor mass and EDEP in Figure 4(a) and for different families of progenitor mass and f in Figure 4(b). In these figures, a shape scheme is used to denote models with different ZAMS mass, while a color scheme is used to denote families of models with constant EDEP and constant f in Figures 4(a) and (b), respectively. In Figure 4(b) we note that models that correspond to a fixed progenitor mass and f fall in a curve that displays an inverse correlation between luminosity and FWHM. Points in a cluster with the same shape and color are characterized purely based on their EDEP. By comparing the points in such clusters with the corresponding points in Figure 4(a), we infer that the higher EDEP models give rise to higher peak luminosities and lower FWHMs, just as in the case of Dutch mass-loss models. From Figure 4(b), we also infer that clusters of models with a larger progenitor mass lie toward the lower end of the FWHM axis with the exception of 13 M e , f = 0.015 models.
In Figure 4(b), when we move across families of models with fixed ZAMS mass (say 13 M e ) and varying f, we observe that models with larger f tend to have a smaller breakout pulse width, unlike the precursor peaks of the Dutch mass-loss models. This trend is heavily influenced by the shock−CSM interaction in the enhanced mass-loss models, due to the presence of a highly dense CSM. Such interactions broaden the post-peak light curve, giving rise to pulse widths that are an order of magnitude larger than those of the Dutch mass-loss precursor peaks.

Comparison of Explosion Characteristics in Optically Thin CSM with Analytic Models
In Section 5.2, we presented analysis involving the FWHM of luminosity peaks. We also presented a discussion on the method adopted for calculating FWHM from a Gaussian fit to the luminosity peaks. In the comparison with the Calzavara & Matzner (2004) analytical expectation for breakout pulse duration (t se ), we find that FWHMs found by fitting Gaussians to the precursor peaks until the first minimum in L bol that occurs after the precursor peak give a better correlation with t se , rather than fitting Gaussians by omitting the intermediate plateau phase as described earlier.
Here we study the validity of Equation (6) of Calzavara & Matzner (2004) for bolometric luminosities and blackbody temperatures at precursor-peak times for either all of the Dutch mass-loss models or a subset of the 14 M e models that obeys the relation peak µ T L t or 1 BB 4 Peak se . The fit to precursor-peak luminosities is displayed in Figure 5. Gaussian-fitted FWHMs and the retarded times in STELLA at which precursor peaks occur are exploited as candidates for t se or time to shock emergence.
In cases where the different families of overshoot parameter offer much less scatter in the plots above, such as in Figure 5(a), a fit of q/t, where q is a fit parameter related to the constant in Equation (6) of Calzavara & Matzner (2004), offers a good match with the observed trend in data points. The fit to q is 8.3 × 10 3 ± 1.2 × 10 2 . In general, when different families of constant overshoot parameter present a dominant scatter in the plots such as in Figure 5(b), an additional parameter "h" is introduced to correctly model this scatter. See Figure 5 for the respective fit functions. Units of the slopes and the additional parameter are the same as the y and x fit variables, respectively.
We performed similar fits to different families of T BB 4 , where each family was characterized by a constant overshoot parameter. We fit for T 4 = d/(t − e), where t is the retarded time in STELLA. The slope "d" varies from 6.96 × 10 4 to 2.12 × 10 5 across the different f families. In appropriate units (discussed below), the spread in "g" and "h" (see Figure 5(b)) is [2.64 × 10 3 , 4.72 × 10 3 ] and [4.54 × 10 2 , 7.43 × 10 2 ], respectively.
The value of the constant in Equation (6) of Calzavara & Matzner (2004) is 7.5 × 10 25 · K 4 · s, invariant of the choice of cgs or SI units, where a is the radiation constant, κ = 0.34 cm 2 g −1 , f (an additional correction factor introduced by Calzavara & Matzner 2004, hereafter referred to as f C−M ) is taken to be 1, and n = 3/2. For the fits to T BB 4 , the expected value of the slope is (after suitable rescaling) 7.5 × 10 5 . This value lies outside the observed range of fitted slopes, 6.96 × 10 4 to 2.12 × 10 5 , and is larger than one standard deviation from the upper bound of the fitted slopes by a factor of ≈2.8. For the fits to peak L bol , we observe the calculated slopes to agree with the observed fits within a factor of 10. This discrepancy may partly be due to the assumption of homologous expansion in the analytic models of Calzavara & Matzner (2004).

High-velocity Ejection of a Low-mass Shell in a Rarefied CSM Established by Dutch Wind
We present the velocity of the different zones of the star and CSM structure (fed by Dutch mass loss from the pre-SN star) as a function of the optical depth ahead of that zone, analogous to Figure 2 of Kozyreva et al. (2020), in Figure 6. While values for velocity (the ordinate) were read directly from mesa.swd, optical depth τ (the abscissa) ahead of a zone at radius r was obtained indirectly from (Rosseland) opacity (κ) and density (ρ) by ò k r dr r R · · , where R is the outermost zone of the star + CSM structure. We employ a color scheme to display profiles at different times, from before the precursor peak to times after the main peak. The intersection of these profiles with the τ = 1.2c/v line gives the properties of the breakout shell, i.e., the shell from which radiation escapes at SBO (Shussman et al. 2016).
The region in these plots that lies between typically 10 −2 and 10 2 of τ and given by colored dashed lines corresponds to the separation between the outermost zone of the star and the innermost zone of the CSM. Since this region is not sampled by any zones in STELLA, we overcome this numerical artifact by first linearly interpolating τ and v in this region, followed by solving for the intersection of τ and 1.2 · c/v. At all of our sampled time instants except the very late 3.0 days, we observe a sharp peak in velocity, between 60 and 76 × 10 3 km s −1 , at very small values of τ, lesser than 10 −2 . The sharp peaks in velocity are those of shells of CSM that are accelerated to high velocities and labeled by different times. The mass contained in this shell is very small (M shell ∼ 5 × 10 −7 M e ). We note that, even at very early times like 0.02 days in the Dutch mass-loss case, L bol has already climbed to one-third of the precursor-peak luminosity (see Figure 2(a)). Such high luminosities ahead of the 1.2c/v surface can produce high acceleration in rarefied gas, which may be the cause of these high-velocity shells. Radiative or shock acceleration of the shells of the CSM depends strongly on their mass densities. For CSMs generated by an additional enhanced mass loss, with densities higher by at least 2 orders of magnitude, velocities achieved by the CSM are much smaller (see Figure 7(b)). As time goes on, the velocity of the shell slows down after it runs into other shells previously ejected from the progenitor star. These high-velocity shells may thus exist only for relatively early times. Their detections would therefore depend on very early-time spectroscopic observations of faint objects against very bright continuum (from the SN itself).
In Figure 6(b), we reproduce the same profile as in Figure 6(a), but with a discontinuous τ-axis to highlight all the significant qualitative features of these profiles in comparison to Figure 2 of Kozyreva et al. (2020). With the exception of 3.0 days, we note that the effect of the different time slices is noticeable only toward the outer regions of lower τ. With an increase in time, the velocity peaks increase and move outward to lower and lower values of τ, as also observed by Kozyreva et al. (2020). We note that the breakout shell is located at τ ≈ 15 except at 3.0 days, at which τ at intersection is ≈47.
We emphasize here that the region that contains the intersection of v and 1.2 · c/τ (at all sampled time instants before 3.0 days) corresponds to the zone boundaries at which the CSM was appended to the star when STELLA was initialized. This is also the region that we refer to as "poorly sampled" by STELLA. The breakout shell was then found to lie toward the outer edge of this poorly sampled region. At 3.0 days, this region lies even inside of the surface that was the star -CSM boundary at the time of initialization of STELLA. This  indicates that at late times, in the spherical phase of expansion, the observed luminosity is dominated by deeper and deeper layers of the star.

Profiles of Radiation Hydrodynamic Variables and Their Time Evolution
We calculated the time evolution of Mach number for a CSM established by enhanced mass loss using the local sound speed. Until 3 days, the outermost zones in the CSM remain subsonic, characterized by a Mach number that is below unity. As can be seen in Figures 7-8, until 3 days, the luminosity shell determined by τ = 1.2 · c/v lies inside the shock. This results in a flat L bol profile until about 3 days in Figure 2(b). Breakout occurs between 3 and 4 days and is characterized by acceleration to supersonic velocities of all shells in the CSM. Photons from the luminosity shell escape ahead of the shock, giving rise to a luminosity flash.

Explosion inside an Optically Thick CSM
We present the early-time evolution of radiation hydrodynamic variables throughout the star and the CSM for the 12 M e , f = 0.020 star with enhanced mass loss prior to collapse and with EDEP = 1.0 FOE in Figures 7-8. The curves labeled "STELLA initialization" in these profiles are obtained from mesa.hyd containing structural information, which, along with mesa.abn, are input files to STELLA. Just before handover to STELLA, the profiles for the CSM are written to mesa.hyd. While the CSM density is explicitly calculated from the mass-loss history, all other variables are initialized with a constant value in the CSM, the constant being the value these variables obtain at the outermost zone of the star. The CSM variables are displayed toward the outermost radial coordinates, from the 271st zone. For all time instants displayed in these profiles, STELLA has completely relaxed to appropriate values of these variables in the CSM, obtained through hydrodynamic modeling. This is checked through explicit early-time analysis of these profiles, soon after initialization of STELLA.
As noted before, Shussman et al. (2016) interpret the τ = 1.2 · c/v surface as the one from which observed luminosity is emitted. From Figure 2 of Kozyreva et al. (2020), we note that τ at this surface is ≈80 at SBO. For our enhanced mass-loss model, τ for this surface is also of the order of 10 2 .
The light curve for this model is given in Figure 2 (b). For this analysis, all our references to L bol and T eff are with respect to Figure 2(b), and our references to variables ρ, v, T, and L are with respect to Figures 7-8. Evolution of other hydrodynamic variables such as optical depth ahead of a zone (τ) and pressure (P) was also studied.
At t = 2 days, which corresponds to the initial flat portion in L bol and T eff , we note that τ ahead of the shock is high, of the order of 10 2 . From Figures 7-8, we infer that the luminosity shell determined by τ = 1.2c/v lies well inside the radial coordinate of the shock determined by examining Mach numbers until 3 days in STELLA time. At 3 days in Figure 7(b), the sharp rise in velocity of the few outermost zones of the CSM is a signature of radiative acceleration. We observe several distinct features at 4 days, which marks the instant just after the peak in L bol , T eff , and T BB . SBO has occurred by this time, which is best described by Figure 7(b). The velocity of the outermost zones of the CSM gets accelerated by the shock to high velocities between 3 and 4 days. The accelerated zones also display a sharp rise in T and L during this interval. At 4 days, τ, ρ, and P fall steeply ahead of the luminosity shell, allowing photons that were trapped behind the shock to escape out of the CSM in a luminosity flash. At 6 days, the breakout flash is almost over, as the majority of the trapped photons have escaped ahead. L of the τ = 1.2c/v shell saturates. The outer layers of the CSM expand outward with We observe several small peaks in T, inside the τ = 1.2 · c/v surface at early times. However, by peak time, at 4 days, these inner peaks become insignificant. These multiple peaks are caused by the interaction of the shock with the highly dense CSM, which increases the thermal energy of matter in the CSM. Significant contributions from shock−matter interactions in the highly dense CSM of the enhanced mass-loss models give rise to a longer or brighter breakout flash, reflected as a broad peak in L bol and T eff .

Explosion inside an Optically Thin CSM Established by Dutch Mass Loss from the Progenitor Star
We now analyze SBO in our Dutch mass-loss, 12 M e , f = 0.025, EDEP = 1.0 FOE model through stellar profiles presented in Figures 9-10. Nakar & Sari (2010) take the timescale for the commencement of spherical expansion as t s = R star /v 0 ∼ 0.17 days. The precursor and main peaks for the Dutch mass-loss model under consideration are separated by 0.082 days, which can be confirmed from a visual inspection of Figure 2(a). Therefore, both precursor (L 1 ) and main (L 2 ) peaks are in the planar phase of expansion.  The orange curves in all of the Dutch mass-loss profiles except L were plotted from the mesa.hyd file, one of the files that STELLA is initialized with. We display L bol and T eff in Figure 2(a), and profiles of variables ρ, v, T, and L at different times are given in Figures 9-10.
The precursor peak in the Dutch mass-loss models is due to SBO from the star. The main peak in L bol and T eff can be explained by the interaction of the forward shock with dense regions of the CSM. This claim can be supported with the following calculation. The peak velocity at precursor-peak time, which is a proxy for shock velocity, for our Dutch massloss model is log(v [cm s −1 ]) ≈ 9.8. The precursor and main peaks are separated by 0.082 days (see Figure 2(a)). The approximate distance that the shock could travel in 0.082 days is about 0.45 × 10 14 cm.
From Figure 9(a), we find the separation between the outer edge of the exploding star and the largest density peaks in the CSM to be of the same order that the shock would travel in the intervening time period between the precursor and the main peaks. In other words, the timing of the main (L 2 ) peak can be rationalized in terms of shock travel time to the nearest dense shell in the CSM.

Comparison of Simulations with Analytic Formulation:
Shock Breakout from an Optically Thick CSM In this subsection, we compare the analytical results of Chevalier & Irwin (2021) in a reduced form (assuming that the opacity-dependent factor on the right-hand side of Equation (3) is unity) with our numerical results obtained via STELLA. This form eliminates the dependence on their parameter D * , which is related to the CSM density profile. We find that the product of the peak luminosity with the square root of the rise time of the luminosity of the breakout obtained from STELLA runs is actually proportional to the predictions of Chevalier & Irwin (2021), i.e., the products of ejecta mass -M ej 1 2 and the outer radial extent R w15 1 2 of the dense CSM established by heavy mass loss from the progenitor star prior to the SN explosion. Here we approximate the rise time by half the value of the FWHM of the luminosity light curves of all relevant cases computed by STELLA for enhanced-mass-loss-fed CSM. Given the nearsymmetric light-curve shape (in time) as seen in Figure 2(b), this is a good approximation. However, to obtain a near one-toone scaling of the analytical and numerical results, we had to use an overall multiplicative factor 0.66 to the reduced relation obtained from Equation (3) (see Figure 11(a)).
In Figure 11(b), we select the cases of two specific energies of explosion (EDEP = 0.8 and 1.2 FOE) among all the runs with different ZAMS masses and overshoot factors f and display the ratios of the analytical products of the peak luminosity L Peak and t rise , multiplied by a factor of 0.66, and the corresponding STELLA results as a function of f. We find that the ratios are relatively close to a factor of unity except in an intermediate range of f values. The convective overshoot factors affect stellar structure and evolution and influence the timescales a star spends in a given part of the Hertzsprung-Russell (H-R) diagram. Since that, in turn, controls the massloss rate from the progenitor star and the ejecta mass left on the star at the time of explosion, and even the farthest extent of the wind-fed dense CSM, it is not surprising that the overshoot factors may have signatures that are not fully captured in the analytical formulation. We note that for some of these progenitor ZAMS masses the range 0.22 f 0.36 for the overshoot parameter is where the stars are susceptible to undergoing blue loops in the H-R diagram as opposed to staying put in the RSG regions until core collapse and explosion. Hence, while the overshoot factors are known to affect the size and structure of the star, they may even influence the characteristics of the dense CSM and therefore also the SBO process, which require more detailed modeling.

Photon Spectral Energy Distribution of Shock Breakout of a Supernova
We illustrate here the temporal evolution of the spectral energy distribution ν · L ν following Blinnikov & Tolstov (2011) for our prototypical Dutch mass-loss model in Figure 12 and for our prototypical enhanced mass-loss model in Figure 13. An interesting feature in our spectra is the occurrence of double peaks in Figure 12, which become prominent at late times, unlike the single-peaked profile in Figure 9 of Shussman et al. (2016). By studying the time evolution of these profiles, we infer that the double peaks originate at a time immediately after the precursor (L 1 ) peak and hence seem to be heavily influenced by SBO itself. This double-peaked spectral profile is present in several progenitors of various masses and other attributes, but only with Dutch mass loss  m. The enhanced mass-loss models display only a single-peaked profile.
We identify the Rayleigh-Jeans regime by the condition h · ν < k B · T obs following Shussman et al. (2016), where we take the blackbody temperature, T BB reported by STELLA in the mesa.tt file, as T obs . Shussman et al. (2016) note that deviations from the blackbody spectrum, L ν ∼ ν 2 , in the Rayleigh-Jeans regime are significant during the breakout and planar phases and suggest the correction L ν ∼ ν 1.4 . We define two functions n n L log( · ) = u + n log 3 · ( ) and n n L log( · ) = w + n 2.4 log · ( ) and fit these functions for optimal values of u and w that give the best fit to spectral data in the Rayleigh-Jeans regime. Here typical values of u and w are −2 and 6.5 for the enhanced mass-loss cases and −3.6 and 5.2 for the Dutch mass-loss cases, respectively. We show these fits for ν in the Rayleigh-Jeans regime for each of the panels in Figures 12 and 13. In panels (b) and (c) of Figure 12 for the spectra of Dutch mass-loss cases, we note a second spectral bump well beyond the frequencies of the first peak. The time evolution of blackbody temperature for this model can be found in Figure 2(a), where T BB reached a relatively large value.
We observe from Figure 12 that the blackbody spectrum fit for our flux distribution in the Rayleigh-Jeans regime better than the ν 1.4 correction at early-time instants, until the mainpeak time, by which we observe a significant peak in T BB in Figure 2(a). The ν 1.4 correction offers a better fit from mainpeak time. We also note that in panel (b), at main-peak time, the blackbody fit is better except for frequencies close to the first spectral peak. At very late times such as 0.7 days in panel (c), the ν 1.4 correction clearly offers a better fit. In Figure 13, we observe noticeable differences in the two fits only from peak time in panel (b). From peak time, in panels (b)-(c), we  observe the ν 1.4 correction to clearly give a better fit to the spectra in the Rayleigh-Jeans regime.

Discussion
In a previous work (Das & Ray 2017), we have modeled progenitor stars with MESA and simulated their explosions with the SuperNova Explosion Code SNEC. 12 SNEC is a spherically symmetric 1D Lagrangian radiation hydrodynamics code that simulates SN explosions, including the propagation of the shock through the stellar envelope, and produces bolometric and multicolor light curves. A limitation of SNEC is the assumption of local thermodynamic equilibrium (LTE), and the same temperature is imposed for radiation and matter. This is problematic especially at the time of SBO and during the transitional stage when the ejecta goes from optically thick to optically thin states (Morozova et al. 2015). SNEC utilizes essentially gray (i.e., frequency-integrated) radiation hydrodynamics, whereas the state-of-the-art codes like CMFGEN (Dessart et al. 2013) solve non-LTE line transfer radiative transfer. STELLA handles multigroup radiation hydrodynamics (Moriya et al. 2011) and implements radiative transfer without assuming an equilibrium photon spectrum (Baklanov et al. 2005).
The special care in utilizing the appropriate opacities by the developers of STELLA is one of its many strong suits over analytic and gray-transport calculations using simplifying but limiting assumptions. Blinnikov & Bartunov (1993) point out that one-group calculations need to distinguish between energy-averaged opacity κ E and the Planck mean 13 κ P in the energy equation and between the flux mean κ F and the Rosseland mean 14 κ R . In the literature, however, many authors have simply put κ E ≡ κ P and κ F ≡ κ R as also κ E ≡ κ F ≡ κ R . The latter assumptions are inadequate to handle situations such as when hot ultraviolet radiation falls on cold static matter (Marshak 1958;Blinnikov & Bartunov 1993). In such cases, the temperature inside a medium may become sufficiently large to generate a radiation wave, 15 which in turn can affect the propagation of the shock into a second medium bounding the first (Marshak 1958). The multigroup radiation hydrodynamics in STELLA is free of this limitation, and the code correctly accounts for the large changes in monochromatic opacity over many orders of magnitude (for example, at the Lyman jump).
The closely spaced precursor and following bolometric luminosity pulse(s) found in our simulations occur very early after SBO. Thus, it may be difficult to detect observationally these flashes unless one is forewarned about the arrival of the shock wave at the surface, such as through the arrival of a neutrino signal (Antonioli et al. 2004) from the core bounce and delayed shock revival in an SN. These events necessarily have to be relatively close to the Earth for effective neutrino detection and therefore have low event rates. In addition, after the first detection, follow-up (multiband) observations also have to be carried out soon afterward. However, the advent of very wide field deep searches for optical transients (e.g., LSST Science Collaboration et al. 2009;Rau et al. 2009;Graham et al. 2019) has made possible the early detection of SNe, e.g., as in SN 2018fif. This arose from an RSG progenitor in UGC 85 at a distance of about 73 Mpc with an r-band detection within 0.2 days after a "reference time" t 0 at which point the flux was zero (Soumagnac et al. 2020). Hence, the era of such early detections is already upon us, and it may be simply a matter of time before detection of peculiar features in very early light curves will be routinely carried out. Nakar & Sari (2010) developed an analytic model based on the Sakurai (1960) model to provide luminosity and temperature evolution ("light curves") starting from breakout, through the planar and spherical expansion phases, until radioactive decay input or recombination in a cooler envelope starts to play a substantial role. They described an example of SBO of an RSG that had an early bright peak in optical and UV bands, less than an hour after breakout, followed by a minimum at the end of the planar phase and then another peak again once the color temperature of the pulse drops to the observed frequency range. This is qualitatively similar to the detection and use of the light flash of an atmospheric nuclear explosion (which has a characteristic double-peaked light pulse) to measure the yield of the explosion. The times of the minimum and second maximum were shown to be related to the square root of the yield (Barasch 1979). The second peak in such atmospheric nuclear tests lasts about 100 times longer than the first peak but contains almost 99% of the total radiated energy, which is about one-fourth of the yield. The rapidly expanding shock front resulting from the explosion is initially opaque to visible light, but as expansion proceeds, the shock front continues to cool and becomes less opaque to visible light, allowing light from the hotter interior to begin to escape and leading to the increase of optical intensity again 16 (Barasch 1979).

Conclusions
In this paper we have investigated the characteristics of SBO radiation with the help of a multigroup radiation hydrodynamics code STELLA starting from a suite of pre-SN models that we have evolved from ZAMS until core collapse with the 12 https://stellarcollapse.org/SNEC.html 13 The Planck mean yields the correct value for the integrated thermal emission from an optically thin plasma. It is calculated from k r n T , ABS ( ) and the Planck function B ν (T), and its weighting function peaks at hν = 2.8kT, which is where the monochromatic opacity is most strongly sampled when taking the Planck mean (Fontes 2015). 14 The Rosseland mean opacity yields the correct value of the integrated energy flux for an optically thick plasma. It uses the weighting function ∂B ν (T)/∂T. Here the weighting function peaks at hν = 3.8kT. Around this peak the monochromatic total opacity k k k = + n n n TOT ABS SCATT will be most strongly sampled when the Rosseland mean is computed (Fontes 2015). While the scattering part might dominate in a large part of the relevant (ρ, T) phase, it is through the absorptive processes that photon numbers and their energies can change. Hence, the latter is crucial to the approach to thermal equilibrium between matter and radiation. 15 In a similar process of SBO out of a neutrinosphere in a core bounce SN a short precursor neutrino wave is present ahead of the shock at high densities, which lengthens as the shock approaches the neutrinosphere at lower densities, and ultimately the neutrinos escape to infinity as the shock crosses the neutrinosphere (Bethe et al. 1980). 16 Note, however, that both sets of double-peaked flashes (the calculations of Nakar & Sari 2010;Shussman et al. 2016, as well as in the atmospheric nuclear explosions) refer to specific radiation bands, say, optical V, g, or R bands or FUV band, whereas our results described in Section 5 refer to bolometric light curves. Development of multiple peaks of matter density as a function of radial or mass depth from the surface at around the time of breakout are also reported -see Figure 2 of Blinnikov & Tolstov (2011). Closely separated double peaks of bolometric luminosity near breakout time in simulations of Population III core-collapse SNe are also reported (Whalen et al. 2013). In "Eye Witness Report of Nuclear Explosion, 16.7.45," Frisch (1964) describes the temporal variations of the shape, brightness, and color of the light source in the explosion and ascribes the sudden appearance at late times of multiple white patches on the underside of the cloud layers just above the explosion to the impact of the blast wave on the cloud layers in the surrounding medium. These appeared before the initially spherical (and red) fireball had started to flatten out. help of the MESA code. Our pre-SN stars have been evolved with a variety of inputs on its internal convection and mixing, as well as two sets of mass-loss schema until the core-collapse stage. Finally, we take these pre-SN stars through the corecollapse stage, initiation of a hydrodynamic shock via a thermal bomb mechanism, and follow the shock propagation throughout the star until the shock nearly reaches the surface of the star with the added features in MESA version r-10398, until finally it is handed over to the STELLA code to model the SBO stage with a multigroup radiative transfer scheme coupled with Lagrangian hydrodynamics of the explosion.
Our main findings in the paper can be summarized as follows. First, we find closely time-separated, double peaks in bolometric light curves computed by the STELLA code and explain these as a combination of SBO flash and the immediately following interaction of the shock wave with shells of dense gas in the CSM ejected earlier during pre-SN evolution.
Second, we find the presence of high-velocity ejecta, containing minuscule mass, that was earlier ejected in the Dutch wind but was accelerated during the SBO and flash. This is a qualitatively different type of ejecta velocity structure than that of the bare star SBO considered by Kozyreva et al. (2020), Shussman et al. (2016), and others, even in the Dutch wind CSM, which is optically thin.
Third, we show how the breakout flash from an optically thick CSM differs from that of the optically thin CSM case fed by pre-SN Dutch wind scenarios. Moreover, we find that the photon spectrum of the breakout flash, although it appears to have an initial thermal shape, would very soon develop into a spectrum that is clearly far from a thermal shape. It has a large, disjoint high-energy bump visible from a separate layer generating these energetic photons. This is seen in the case of the optically thin CSM fed by the Dutch wind scheme and is absent in the photon spectrum emerging from an optically thick CSM.
We also place our results of numerical simulation of the SBO properties in the context of analytical results such as those by Calzavara & Matzner (2004) for thin wind and Chevalier & Irwin (2011, 2021 for optically thick wind. We compare the characteristics of the ensemble of SBO flashes out of our suite of pre-SN models with the analytical predictions of peak bolometric luminosities and rise time of the breakout flash under some simplifying approximations.