Gamma-ray Thermalization and Leakage from Millisecond Magnetar Nebulae: Towards a Self-Consistent Model for Superluminous Supernovae
DDraft version January 15, 2021
Typeset using L A TEX twocolumn style in AASTeX63
Gamma-ray Thermalization and Leakage from Millisecond Magnetar Nebulae:Towards a Self-Consistent Model for Superluminous Supernovae
Indrek Vurm and Brian D. Metzger
2, 3 Tartu Observatory, Tartu University, 61602 T ¯ o ravere, Tartumaa, Estonia Department of Physics and Columbia Astrophysics Laboratory, Columbia University, Pupin Hall, New York, NY 10027, USA Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave, New York, NY 10010, USA
ABSTRACTSuperluminous supernovae (SLSNe) are massive star explosions that are too luminous to be poweredby traditional energy sources, such as the radioactive decay of Ni. These transients may instead bepowered by a central engine, such as a millisecond pulsar or magnetar, whose relativistic wind inflatesa nebula of high energy particles and radiation behind the expanding supernova ejecta. We presentthree-dimensional Monte Carlo radiative transfer calculations of SLSNe which follow the productionand thermalization of high energy radiation from the nebula into optical radiation and, conversely,determine the gamma-ray emission that escapes the ejecta without thermalizing. Specifically, we trackthe evolution of photons and matter in a coupled two-zone (“wind/nebula” and “ejecta”) model, ac-counting for the range of radiative processes including (typically multiple generations of) pair creation.We identify a novel mechanism by which γγ pair creation in the upstream pulsar wind regulates themean energy of particles injected into the nebula over the first several years after the explosion, ren-dering our results on this timescale insensitive to the (uncertain) intrinsic wind pair multiplicity. Toexplain the observed late-time steepening of SLSNe optical light curves as being the result of gamma-ray leakage, we find that the nebular magnetization must be very low ε B ∼ < − − − . For higher ε B synchrotron emission quickly comes to dominate the thermalized nebula radiation, and being readilyphotoelectrically absorbed because of its lower photon energies, results in the supernova optical lightcurve tracking the spin-down power even to late times ∼ > ∝ t − in a way that coincidentally mimicksgamma-ray escape. Our models predict a late-time flattening in the optical light curves of SLSNe, forwhich there may be evidence in SN 2015bn. Keywords: supernovae, gamma-rays INTRODUCTIONA growing sample of stellar explosions have been dis-covered with luminosities too high to be powered bytraditional energy sources, such as the radioactive de-cay of Ni. These include the rare class of stellar corecollapse events known as “superluminous supernovae”(SLSNe; e.g. Quimby et al. 2011; Nicholl et al. 2013; In-serra et al. 2013; Howell et al. 2013; Dong et al. 2016; DeCia et al. 2018; Lunnan et al. 2018; Quimby et al. 2018;see Gal-Yam 2018; Inserra 2019 for recent reviews) aswell as the emerging class of luminous transients withfaster evolving light curves indicative of a lower ejectamass, sometimes referred to as “fast blue optical tran-sients” (FBOTs; e.g. Drout et al. 2014; Arcavi et al. 2016; Tanaka et al. 2016). The best studied FBOT isAT2018cow (Prentice et al. 2018; Perley et al. 2019;Kuin et al. 2019), which peaked on a timescale of onlya few days and was accompanied by nonthermal emis-sion across the electromagnetic spectrum, from radio togamma-rays (Margutti et al. 2019; Ho et al. 2019).Powering the light curves of SLSNe and FBOTs re-quires prolonged heating of the ejecta by a centrallyconcentrated energy source. One such potential sourceare shocks generated by the collision of the ejecta witha dense circumstellar shell or disk surrounding the pro-genitor star at the time of explosion (e.g. Smith & Mc-Cray 2007; Chevalier & Irwin 2011; Moriya et al. 2013).Circumstellar interaction is most compelling as an ex-planation for the hydrogen-rich class of SLSNe (which a r X i v : . [ a s t r o - ph . H E ] J a n Vurm & Metzger often show emission lines indicative of slow-moving gas;e.g. Smith et al. 2007; Fransson et al. 2014; Benetti et al.2014; Nicholl et al. 2020), though similar shock-poweredinteraction could also account for the light curves of atleast some hydrogen-poor SLSNe-I (e.g. Sorokina et al.2016; Kozyreva et al. 2017; Rest et al. 2018).Another commonly discussed model for SLSNe-I/FBOTs invokes energy injection by a young activecompact object, such as an accreting black hole or neu-tron star (Quataert & Kasen 2012; Woosley & Heger2012; Margalit & Metzger 2016; Moriya et al. 2018), orthe rotationally-powered wind of a pulsar or magnetarwith a millisecond spin period (Maeda et al. 2007; Kasen& Bildsten 2010; Woosley 2010; Dessart et al. 2012; Met-zger et al. 2015; Sukhbold & Woosley 2016; Sukhbold &Thompson 2017). Among the evidence supporting thepresence of a central engine are emission line features inthe nebular spectra of SLSNe (e.g., Nicholl et al. 2019)and peculiar Type Ib SNe (e.g., Milisavljevic et al. 2018)similar to those seen in the hyper-energetic supernovaeobserved in coincidence with long duration gamma-raybursts. A central engine in AT2018cow is supported bythe rapid time variability and energy spectrum (particu-larly the Compton hump feature) of the X-rays observedin coincidence with the optical radiation (e.g. Marguttiet al. 2019).Theoretical models for the visual light curves andspectra of engine-powered supernovae generally assumethat the power output of the central engine is ther-malized by the surrounding ejecta with 100% efficiency(Kasen & Bildsten 2010; Woosley 2010; Dessart et al.2012; Chatzopoulos et al. 2013; Metzger et al. 2015).While this assumption is important to the overall vi-ability of the scenario, and to quantitative inferencesdrawn about the properties of the central engine (e.g.,the dipole magnetic field strength and birth spin periodof the neutron star), its justification has thus far receivedlittle attention in the literature.A broad outline for how such a thermalization processmight occur is as follows. A relativistic wind from theneutron star or black hole engine inflates a nebula ofrelativistic electron/positron pairs and radiation behindthe expanding supernova ejecta shell (Kotera et al. 2013;Metzger et al. 2014; Murase et al. 2015). Pairs heatednear the wind termination shock enter the nebula andquickly radiate their energy via synchrotron and inverseCompton (IC) processes in a broad spectrum spanningthe X-ray/gamma-ray band, in analogy with ordinarypulsar wind nebulae like the Crab Nebula (B¨uhler &Blandford 2014). Only the portion of this radiationwhich can escape the nebula without undergoing γγ paircreation, and is in the appropriate energy range to be absorbed and thermalized by the ejecta shell, is availableto heat the ejecta and power the supernova emission.Thermalization of the nebular radiation will be mostefficient at early times after the explosion, when the col-umn density through the ejecta shell and “compactness”of the system are at their highest. However, as theejecta shell expands and becomes increasingly transpar-ent its radiation field dilutes. As a result, the efficiencyof the thermalization process will drop and the opticalsupernova light curve will eventually fall below the rateof energy injection from the central engine, as the ma-jority of the engine’s power escapes directly from thenebula as high-energy radiation.In support of such a picture, the well-studied SLSNe2015bn (Nicholl et al. 2016a,b) showed a marked steep-ening of its optical/UV light curve at t ∼ >
200 days be-low the predicted ∝ t − decay of the magnetar dipolespin-down luminosity, qualitatively consistent with theexpectation of radiation leakage (Nicholl et al. 2018).Searches for this “missing energy” by means of X-ray(Margutti et al. 2018; Bhirombhakdi et al. 2018) andgamma-ray (Renault-Tinacci et al. 2018) observationsof SLSNe at late times have thus far only resulted in up-per limits (with one possible exception; see Levan et al.2013). However, this is not necessarily surprising be-cause the ejecta is expected to remain optically thickto soft X-rays for decades after the explosion (Margalitet al. 2018a), while most of the existing gamma-ray up-per limits are not deep enough to constrain the expectedescaping flux.In this paper, we present calculations of engine-powered supernova light curves which for the first timeaccount self-consistently for the thermalization and es-cape of high energy radiation from the central nebulaand the surrounding ejecta shell. We accomplish thisby means of three-dimensional time-dependent MonteCarlo simulations which track the coupled evolution ofphotons and electron/positron pairs in the nebula andejecta through the myriad physical processes couplingthem. Such a detailed treatment is necessary becausethe process of gamma-ray thermalization is complex andnon-linear.For example, although many processes in the su-pernova ejecta shell and its radiation field can absorbgamma-ray photons, the initial absorption of a photonis no guarantee of its ultimate thermalization. Sev-eral photon destruction processes result in the creation At early times, the nebula can also deposit energy directly bydriving a shock into the inner edge of the ejecta shell (e.g. Metzgeret al. 2014; Kasen et al. 2016; Chen et al. 2016; Suzuki & Maeda2019, 2020). amma-ray thermalization and leakage from millisecond magnetar nebulae γγ pair creation also occurs in the up-stream region of the pulsar wind (interior to the termi-nation shock), augmenting the wind mass-loading andregulating the average energy of the pairs the wind re-leases into the nebula. Thus, even the inner boundarycondition defined by the engine is intimately coupled tothe properties of the larger-scale nebula/ejecta duringthe most accessible observational window.This paper is organized as follows. We begin in Sec-tion 2 with several preliminary considerations which mo-tivate the details of our model. These include a briefoverview of engine-powered supernovae (Section 2.1);properties of the radiation field in the nebula (Section2.2); an enumeration of key photon absorption processes(Section 2.3); the process of thermalization of high-energy photons by the ejecta (Section 2.4); and regula-tion of the wind mass-loading by in-situ γγ pair produc-tion (Section 2.5). In Section 3 we describe and presentthe results of our numerical radiative transfer calcula-tions. Further interpretation and implications of ourresults are discussed in Section 4. We summarize ourfindings and conclude in Section 5. The Appendices pro-vide several auxiliary calculations which support thosein the main text. PRELIMINARY CONSIDERATIONSThis section introduces the key concepts and physi-cal processes that enter our calculations and providesseveral analytic estimates that are useful later in inter-preting our numerical results. The reader uninterestedin these details on the first pass is encouraged to advancedirectly to the description and results of the simulations(Section 3), returning to this section as needed. Table 1summarizes several of the key timescales in the problem.2.1.
Overview of engine-powered supernovae
The result of the supernova explosion is the creationof a ballistically expanding shell of ejecta of mass M ej = 10 M M (cid:12) , mean velocity v ej = 10 v cm s − ,mean radius R ej = v ej t , and thickness ≈ R ej . The mean The high pressure of the nebula may drive a shock into the ejectashell and lead to its acceleration at early times (e.g. Suzuki &Maeda 2016). However, this acceleration will subside once radia-tion begins to escape from the nebula at times t ∼ > t pk (eq. 7) ofgreatest interest in this paper. The ejecta kinetic energy definedhere is thus the sum of the initial energy of the explosion andthat released by the engine at t ∼ < t pk (e.g., Metzger et al. 2015;Soker 2016). electron density in the shell, and its radial Thomson op-tical depth at time t = 1 t yr yr following the explosionare given, respectively, by n e (cid:39) M ej Y e π ( v ej t ) m p ≈ . × M v − t − cm − , (1) τ T (cid:39) n e σ T R ej = 3 M ej κ es π ( v ej t ) ≈ . M v − t − , (2)where σ T is the Thomson cross section, κ es = Y e σ T /m p ,and we have assumed an electron fraction Y e = 0 . τ T = 1) after atime t T (cid:39) . M / v − yr . (3)A central compact object is assumed to deposit energyinto a hot nebula behind the ejecta at a rate L e ( t ) = E e t e ( α − t/t e ) α ≈ t (cid:29) t e E e t e ( α − t/t e ) α , (4)where E e is the total energy of the engine and t e the du-ration of its peak activity. The power-law decay index α = 2 is the case of an isolated neutron star or magnetarwith a fixed dipole magnetic field and spin inclination(Spitkovsky 2006) and without significant gravitationalwave losses, α (cid:39) .
38 for a fall-back accreting neutronstar (Metzger et al. 2018b), α = 5 / α ≤ / α = 2). Following the conven-tion used by Nicholl et al. (2017b), the engine energy E e (cid:39) . × ( P / − erg is the magnetar’s initialrotational energy (for assumed mass 1.4 M (cid:12) and birthspin period P ), while the engine lifetime, t e (cid:39) . × (cid:18) P (cid:19) B − s , (5)is the magnetic dipole spin-down timescale, where B d = B × G. The spin-down power at times t (cid:29) t e isthus given by (eq. 4) L e = E e t e /t ≈ . × B − t − erg s − , (6) Vurm & Metzger which is notably independent of P . Fitting the mag-netar model to SLSNe light curve data yield values P ∼ − B ∼ . − t e with the photon diffusion timescale on which thesupernova optical light curve peaks (Arnett 1982), t pk ≈ . M / v − / (cid:18) κ opt κ es / (cid:19) / yr . (7)The latter is set by the condition τ opt ∼ c/v ej , where τ opt = ( κ opt /κ es ) τ T is the optical depth of the ejectato thermal optical/UV photons and κ opt is the effectivecontinuum opacity of Doppler-broadened atomic lines(e.g. Pinto & Eastman 2000).The energy deposition rate L e can also be character-ized by the dimensionless “compactness” parameter (e.g.Guilbert et al. 1983), which we define as (cid:96) inj ≡ σ T L e πm e c R ej ≈ t (cid:29) t e . × − B − v − t − (8)An analogous quantity for the radiation field permeatingthe ejecta can be defined as (cid:96) ≡ u γ σ T R ej m e c , (9)where u γ is the energy density of the radiation field inquestion . The above quantities are related such thatif all the injected energy is converted to radiation inan optically thin environment (with a photon residencetime ∼ R ej /c ), then one would have (cid:96) ≈ (cid:96) inj . The mainutility of the compactness parameter lies in its directrelation to important quantities within the ejecta andnebula such as the pair-production opacity and electroncooling rate, allowing for a readily scalable descriptionof the physical system.A portion of the injected luminosity L e will be ab-sorbed by the ejecta and reprocessed to thermal opti-cal/UV wavelength radiation of characteristic thermal Much of the inferred parameter space includes dipole magneticfields B d ∼ < G outside the range of most Galactic magne-tars, instead overlapping those of radio pulsars. Furthermore,while the usual hallmark of Galactic magnetars is flaring activitydriven by the dissipation of magnetic energy, here we are inter-ested in the comparatively steady extraction of rotational energy,again more akin to pulsars than magnetars. Nevertheless, we usethe term “magnetar” throughout this manuscript to remain con-sistent with the SLSNe literature. Compactness (cid:96) of a radiation field of energy density u γ withina region of size R is numerically equal to the Thomson opacityof the same region that would be obtained if all the radiationenergy was converted to electron/positron rest mass. temperature T eff ∼ − K, which is responsible forpowering the observed supernova emission. It is thususeful to introduce a “thermal compactness” parameter, (cid:96) th ≡ u th σ T R ej m e c = (cid:96) inj f th (1 + τ opt ) , (10)where u th ≈ L e (1 + τ opt ) f th / (4 πcR ) is the energy den-sity of the thermal optical/UV supernova radiation ofcharacteristic photon energy (cid:15) opt ≈ kT eff ∼ −
10 eV,and the factor (1 + τ opt ) accounts for the additional res-idence time of thermal photons in the ejecta.The factor f th entering Equation 10 is the efficiencywith which the engine luminosity is thermalized by theejecta. A self-consistent determination of f th is non-trivial and will take up the greater part of this pa-per. However, for purposes of an estimate connect-ing to the prior literature, we can parameterize f th =1 − exp( − τ therm ) ≈ τ therm via a gamma-ray “thermal-ization” optical depth τ therm ≡ ( κ th /κ es ) τ T (Wang et al.2015; Chen et al. 2015), which does a good job ex-plaining the late-time decay rate of SLSNe light curvesfor values of the effective thermalization opacity κ th ∼ . − . g − (e.g., Nicholl et al. 2017c). Parame-terized thus, one finds that at times t (cid:29) t pk ( f th (cid:28) (cid:96) th ≈ . × − (1 + τ opt ) M v − B − (cid:18) κ th /κ es . (cid:19) t − . (11)One can also define a “nonthermal” compactness param-eter (cid:96) nth of the nebula, which is analogous to (cid:96) th but with u γ replaced by the density of nonthermal photons and R ej with the nebula radius R n .As we shall discuss in Section 2.3, the optical depthacross the nebula and ejecta to γγ pair production tophotons near the threshold energy ( xθ ≈
1) can beroughly expressed as τ γγ, th ∼ (cid:96) th / (15 θ ) ∼ (10 − ) (cid:96) th ,where x ≡ hν/m e c , θ ≡ kT eff /m e c ∼ − − − and T eff ∼ − K is the effective temperature of the tar-get radiation field. The ejecta thus becomes transparent( τ γγ, th ∼ <
1) to photons of energy m e c /θ ∼ − eV on the timescale t γγ ≈ . M / v − / B − / (cid:18) κ th /κ es . (cid:19) / (cid:18) T eff K (cid:19) − / yr , (12)where at these late times we take τ opt <
1. Photons of ∼ TeV energy are thus free to escape to a distant ob-server starting years following the explosion.To connect the above discussion to observations, thetop panel of Figure 2 shows the thermalized luminosity(essentially equal to the optical luminosity after peak),and the inferred escaping engine luminosity, from a sam- amma-ray thermalization and leakage from millisecond magnetar nebulae Figure 1.
Schematic illustration of the engine-powered su-pernova model discussed in this paper. The pulsar wind ofluminosity L e terminates at a shock or region of reconnec-tion (radius R t ), inflating a nebula of relativistic pairs (radius R n (cid:39) R t ) which radiate gamma-rays generated by IC scatter-ing and synchrotron emission. A fraction of the gamma-rayenergy is absorbed and thermalized by the supernova ejecta(of mass M ej , mean velocity v ej , and radius R ej (cid:39) v ej t ) andultimately reprocessed into optical/UV light, while the re-mainder escapes directly to the outside observer. The ther-malized fraction of the spin-down luminosity decreases intime as the column through the ejecta shell (Thomson opti-cal dept τ T ∝ t − ) and the background thermal and nonther-mal radiation fields (e.g. thermal compactness (cid:96) th ∝ t − )decrease. Within the first few years of the explosion, γγ interactions between gamma-rays and optical/X-ray radia-tion load the pre-shock pulsar wind (at radii ∼ < R t ) withelectron/positrons pairs, regulating the flux ˙ N ± and meanenergy L e / ˙ N ± m e c of the pairs which enter the nebula andgenerate the gamma-rays. ple of 38 SLSNe from Nicholl et al. (2017b), who fit ob-served multi-band optical light curve data to a modifiedversion of the Kasen & Bildsten (2010) magnetar modelwithin a Bayesian framework to determine properties ofthe engine ( E e , t e or equivalently B d , P ) and the super-nova ejecta ( M ej , v ej , κ th ). The bottom panel of Figure 2shows the time evolution of the Thomson optical depth, τ T (eq. 2), and the approximate γγ optical depth of a ∼ TeV photon on the supernova optical radiation, τ γγ, th (eq. 30). 2.2. Nebular radiation
This section describes the sources of pairs and radi-ation in the nebula. We again focus on the case ofa millisecond magnetar, but many of our conclusionswould apply to other central engines (e.g. an accretion-
Figure 2.
Top Panel:
Thermal luminosity L th ≡ L e (1 − e − τ therm ) (brown lines) and inferred “escaping” gamma-rayluminosity L esc ≡ L e e − τ therm (black lines), calculated usingejecta and magnetars parameters inferred from Bayesian fitsof optical light curve data to the magnetar model for a sampleof 38 SLSNe from Nicholl et al. (2017b). Here τ therm is theoptical depth of “thermalization”, which we have calculatedassuming a fixed thermalization opacity κ th = 0 .
01 cm g − motivated by a phenomenological light curve model fit toSN2015bn. The specific cases of SN2015bn and SN2017egm,whose properties motivate the numerical models presentedin this paper (Table 3), are shown as red and blue lines,respectively. Bottom Panel:
For the same models in the toppanel, the Thomson optical depth τ T (eq. 2; brown lines)and an estimate of the optical depth of ∼ TeV photons to γγ pair production on the optical supernova light, τ γγ, th (eq. 30;black lines). powered jet) provided that the mean energy per particlefrom their relativistic outflows is similarly high to a pul-sar wind.The characteristic rate of particle injection from a ro-tationally powered pulsar wind is given by (Goldreich & Vurm & Metzger
Table 1.
Summary of Key TimescalesSymbol Description of Event Typical Value † t e (eq. 5) Engine lifetime (e.g., magnetar dipole spin-down time) ∼ . − . t pk (eq. 7) Peak of the supernova optical light curve ∼ . t B (eq. 17) Synchrotron dominates IC cooling in nebula ( (cid:96) B ≈ . (cid:96) th ) ∼ . ε − / , − yr t T (eq. 3) Ejecta transparent to Thomson scattering ( τ T ≈ ∼ t γγ (eq. 12) Wind/nebula and ejecta transparent to γγ pair creation ( τ γγ, th ≈ ∼ t ± (Table 2) Pair-loading of pulsar wind no longer regulated by in situ γγ interactions ∼ −
10 yr (Table 2) t c (eq. 23) Nebula transitions from radiative to adiabatic evolution ∼ (10 − ε / , − yr † For typical values of the ejecta mass M ej ∼ M (cid:12) and mean velocity v ej ∼ − . Julian 1969)˙ N ± ≡ µ ± (cid:18) Ie (cid:19) ∼ × (cid:16) µ ± (cid:17) B P − s − ≈ t (cid:29) t e . × (cid:16) µ ± (cid:17) B − t − s − , (13)where I ≡ πcR η GJ (cid:12)(cid:12) R L , η GJ ≈ Ω B/ πc is theGoldreich-Julian charge density evaluated at the lightcylinder radius R L = P c/ π , P ms = P/ µ ± is the pair multiplicity of the wind (typically, µ ± ∼ < ; e.g. Timokhin & Harding 2019; however, see Be-loborodov 2019, who find a different expression for ˙ N ± in the case of active magnetars).The injected pairs inflate the nebula of radius R n ≈ v n t ∼ < R ej , where the expansion velocity of the nebulais estimated as v n ≈ v ej / E e , v ej , M ej characteristic of SLSNe (Margalit et al. 2019).We can estimate the strength of the magnetic field inthe nebula B n by assuming that the magnetic energy( B / π ) V n , where V n ≈ πR / ε B = 0 . ε B , − of the injected magnetarenergy ∼ L sd t over the expansion time ∼ t , B n ≈ (cid:18) ε B L e tR (cid:19) / ≈ . ε / , − B − v − / t − G . (14)Magnetization values of ε B ∼ − − . u γ with u B = B / π and R ej with R n to define a “magneticcompactness” of the nebula according to (cid:96) B ≡ σ T m e c B π R n ≈ . × − ε B , − B − v − t − . (15)The ratio of the magnetic and thermal compactnesses(eq. 10) is (cid:96) B (cid:96) th ≈ τ opt ) ε B , − M − v (cid:18) κ th /κ es . (cid:19) − t . (16) Once (cid:96) B ∼ > (0 . − (cid:96) th , synchrotron emission will over-take IC scattering as the dominant cooling mechanismfor relativistic pairs in the nebula, with implications forthe photon energy spectrum of the latter. This criticaltransition occurs on a timescale t B ≈ . (cid:18) (cid:96) B (cid:96) th (cid:19) / M / ε / B, − v / (cid:18) κ th /κ es . (cid:19) / yr , (17)which is less than a month after the explosion for ε B ∼ − . For lower ε B (cid:28) − (which, as we shall find,may be needed to explain late-time SLSNe light curves;Section 4.1) the (cid:96) B ∼ (cid:96) th transition occurs only afterseveral years.Pairs enter the nebula at the wind termination shock,where in normal pulsar wind nebulae they are heated.This heating occurs either at the shock front itself (Ken-nel & Coroniti 1984) or by magnetic reconnection inthe striped pulsar wind ahead of the shock (Sironi &Spitkovsky 2009, 2011; Zrake & Arons 2016). Naively,one would expect the freshly injected pairs to acquire amean random Lorentz factor,∆ γ = γ in (cid:39) L e ˙ N ± m e c ≈ . × (cid:16) µ ± (cid:17) − B − t − , (18)where Eq. (13) is used for the pair injection rate ˙ N ± .However, at early times ∆ γ will generally be muchlower than this estimate because of two effects. Firstly,the upstream wind prior to the termination shock isloaded by secondary pairs generated by γγ processes,increasing the effective value of ˙ N ± sharing the pulsar’sluminosity. As we will show in Section 2.5, at early timesafter the explosion this additional pair-loading regulatesthe mean Lorentz factor the pairs entering the nebulato a value γ ∼ − (cid:28) γ in (see also Table 2). Thisenhanced pair-loading eventually ceases once the com-pactness of the nebula drops sufficiently, typically on atimescale of several years to a decade.Another effect can in principle also limit the energyof pairs entering the nebula at late times, though it is amma-ray thermalization and leakage from millisecond magnetar nebulae E < B , then ∆ γ is also limitedby synchrotron cooling during the particle energizationprocess, to a value (e.g. Cerutti et al. 2014) γ rad ≈ (cid:18) e B w r (cid:19) / ≈ . × (cid:18) R t R n (cid:19) / σ − / v / B / t yr , (19)where r e = 2 . × − cm and B w = (cid:20) L e σ ( σ + 1) cR (cid:21) / (20)is the magnetic field strength near the termination of thewind at radius R t , where σ is the wind magnetization(the second line of eq. 19 assumes σ (cid:28) γ ( t ) (cid:39) min[ γ in , γ rad , γ ] , (21)where ∆ γ (cid:39) γ ∼ − at early times while the windexperiences significant γγ pair-loading, before increasingto ∆ γ ∼ γ in , γ rad ∼ < after several years.Upon entering the nebula with energy ∆ γ , the pairscool via synchrotron radiation and IC scattering onlower energy target radiation. In the Thomson regime,IC cooling will dominate synchrotron cooling for as longas the target photon energy density exceeds that ofthe magnetic field (e.g., (cid:96) th ∼ > (cid:96) B for a thermal targetfield). At late times t (cid:29) t B , synchroton dominates ICscattering as the dominant source of pair cooling. Atsome point, however, even synchrotron radiation can nolonger efficiently cool the nebula. The synchrotron cool-ing time is shorter than the nebula expansion time forpairs above a Lorentz factor, γ c = 6 πm e cσ T B t ≈ ε − , − B v t . (22)Equating this to ∆ γ , we see that the nebula will remainfast-cooling until a time t c ≈ (cid:18) ∆ γ (cid:19) / v − B − / ε / B, − yr . (23)At early times t (cid:28) t c the injected pairs radiate all theenergy they receive instead of transferring it to the ki-netic energy of the ejecta via adiabatic expansion. Thenebula’s luminosity will therefore match that of the cen- tral engine for years to decades after the explosion, i.e.on the timescales of relevance to follow-up observations. A pair of Lorentz factor γ ± can produce IC radiationup to the same scale as its own energy, hν IC , max ≈ γ ± m e c ∼ × (cid:16) γ ± (cid:17) eV . (24)By comparison, the characteristic frequency of syn-chrotron radiation from the same pair in the nebulamagnetic field (eq. 14) is significantly lower, hν syn = min[ hν syn , , hν syn , max ] , (25)where hν syn , = hν B γ ± ≈ (cid:16) γ ± (cid:17) ε / , − B − v − / t − keV , (26) ν B = eB n / πm e c is the Larmor frequency, and hν syn , max ≈
160 MeV is the maximum synchrotron fre-quency in the radiation reaction limited case (Guilbertet al. 1983; Cerutti et al. 2014). For the regulated ener-gies of the injected pairs ∆ γ = γ ± ∼ γ ∼ < expectedfor several years after the explosion, we have hν syn ∼ < t ∼ > t B ; eq. 17).2.3. Absorption and energy loss of high-energy photons
Energetic photons deposited into the nebula interactwith both matter and radiation fields as they diffuse out-ward through the nebula and ejecta. The relevant pro-cesses depend on the photon energy and involve bothscattering and absorption, which can lead to the gen-eration of secondary electron-positron pairs and theirradiation. The net effect of these processes is to repro-cess the primary photon energy towards lower frequen-cies. A fraction of this energy eventually reaches thethermal pool and emerges as optical radiation; quanti-fying this thermalization efficiency or optical depth self-consistently is one of the goals of our Monte Carlo sim-ulations. This section summarizes the main radiativemechanisms involved.The main source of opacity for GeV − TeV photons ispair production on soft radiation fields and nuclei in An exception occurs if the nebula magnetization is extremelylow ( ε B ∼ < − ) in which case pairs cool exclusively through ICscattering. Adiabatic losses then set in on a timescale as shortas months: once the IC gamma-rays start to leak out from thenebula the resulting drop in thermalization reduces the targetoptical radiation background, reducing the IC cooling rate in arunaway process. Vurm & Metzger the ejecta (Bethe-Heitler process). The optical depththrough the ejecta of the former is given by (Zdziarski& Svensson 1989) τ ph − mat = 212 π α fs (cid:20) ln(2 x ) − (cid:21) τ T , (27)where α fs (cid:39) / x ≡ hν/m e c (cid:29)
1, and we haveassumed oxygen-dominated composition.Depending on energy, γγ pair production can takeplace on either thermal (e.g., optical/UV) or nonther-mal (e.g., X-ray) radiation fields. For a broad targetspectrum of radiation energy density u γ , a simple ap-proximation for the γγ pair production opacity is givenby (Svensson 1987): τ γγ, nth ≈ η ( α ) x (cid:96) ν (cid:18) x (cid:19) , (28)where η ( α ) ≡ (7 / − α ) − (1 − α ) − / and α is definedso that du γ /d ln x ∝ x α +1 . We have also defined thedifferential compactness of the target radiation field as (cid:96) ν = σ T m e c d u γ d ln x R ej , (29)which is a more precise, frequency-dependent versionof the compactness introduced in eq. (9). The aboveapproximation for τ γγ relies on the availability of targetphotons near the threshold energy 1 /x and is accurateto within 0 .
3% when − < α < θ ≡ kT eff /m e c and compactness (cid:96) th (eq. 10) a simple em-pirical formula for the pair-production opacity can bewritten as τ γγ, th = (cid:96) th x ln(1 + xθ )2 xθ exp (cid:18) − . xθ (cid:19) , (30)which has an error of <
15% at xθ ∼ > .
1. Here x ≈ . θ is the average thermal photon energy in m e c units.Photons can also lose energy by Compton downscat-tering off electrons in the ejecta. An effective opticaldepth for a photon to lose most of its initial energy canbe defined as τ C , eff ≡ t esc t C ≈ ˙ x C x R ej c (1 + τ KN ) ≈ τ T ln(1 + x )1 + 3 x (cid:20) τ T x ln (cid:18) x (cid:19)(cid:21) , (31)where the term (1 + τ KN ) accounts for the enhancedresidence (diffusion) time of the photon in the opaquemedium. Here we have approximated the particle en-ergy loss rate ˙ x C /x ≈ cσ T n e ln(1 + x ) / (1 + 3 x ) whichis accurate to within <
15% at x (cid:29) kT /m e c (where T ∼ − K is the gas temperature), and τ KN ≈ τ T ln(1 + 8 x/ / x (accurate to within 10 −
15% andasymptotically approaching the exact τ KN at x (cid:28) x (cid:29) E ∼ > . τ bf ≈ . × E − τ T , (32)where E keV ≡ hν/ τ T ∼ > ∼ <
10 keV are absorbed for several years afterthe explosion (Margalit et al. 2018a).Figure 3 shows the effective optical depth as a func-tion of the primary photon energy for different Thom-son columns, τ T . From low to high photon energies,the dominant processes responsible for photon energyloss are: bound-free absorption ( E ∼ <
30 keV); Comptondown-scattering ( E ∼ <
100 MeV); Bethe-Heitler photon-matter pair production ( E ∼ <
10 GeV); pair creation onnonthermal radiation (if present); pair creation on ther-mal supernova radiation ( E ∼ > −
100 GeV).Dashed and dotted lines in Figure 3 show differentapproximate analytical results, which may be comparedwith exact results obtained from Monte Carlo simula-tions of photon diffusion (solid lines). In the Compton-dominated range (0 . ∼ < E ∼ <
100 MeV), the simple ana-lytic expression for τ C , eff (eq. 31) can significantly over-estimate the true effective optical depth when τ T > τ eff as a function of photon energy at several snap-shots in time after the explosion, for ejecta and mag-netar parameters fit to the optical light curve data onSN2015bn and SN2017egm (see Table 3). The lowest op-tical depth occurs in an energy window surrounding ∼ τ eff now as a function of time since explosion, fordifferent photon energies. The ejecta becomes transpar-ent ( τ eff <
1) at energies 1 −
10 GeV accessible to
Fermi -LAT by t ∼
100 days, while transparency is delayed foryears at ∼ TeV energies because of γγ absorption by thesupernova optical light.2.4. Gamma-ray thermalization amma-ray thermalization and leakage from millisecond magnetar nebulae E (MeV) e ff T = 30 T = 10 T = 1 T = 0.1MC simulationSpherical diffusion C,eff + abs Figure 3.
Effective optical depth, τ eff , as a function ofinitial photon energy for diffusion out of a spherical cloud,defined such that exp ( − τ eff ) is the fraction of energy thatescapes without significant downgrading (i.e. at photon en-ergies E out > E/ τ T = 30 , , , .
1. Red solid lines show the resultsof full Monte Carlo simulations; dashed lines correspondto an analytical solution to the radiative diffusion equation(Appendix A); and dotted lines show τ C , eff + τ abs , where τ abs = τ ph − mat + τ γγ, th + τ bf is the sum of absorption opti-cal depths and τ C , eff is an estimate of the effective opticaldepth due to Compton down-scattering (eq. 31). The ther-mal compactness in the Monte Carlo calculations is taken tobe (cid:96) th = 10 − , however our results only depend on this as-sumption at the highest photon energies E ∼ > eV (whichundergo γγ interactions on the thermal radiation field). The efficiency of reprocessing and thermalization ofthe nebula energy by the supernova ejecta is crucial toexplaining the optical light curves of engine-powered su-pernovae. However, the mechanism of reprocessing iscomplex and non-linear. Even if gamma-rays are at-tenuated during their escape from the ejecta (Section2.3), this does not guarantee that their energy will bethermalized. In particular, the result of γγ and Bethe-Heitler photon-matter processes (which dominate the at-tenuation of high-energy gamma-rays) is the productionof an energetic electron/positron pair. As shown in Ap-pendix B, this pair will generally experience rapid radia-tive cooling, producing a second generation of radiationwhich itself may or may not escape the ejecta (poten-tially producing yet another generation of pairs, and soon).Energy is deposited into the thermal radiation fieldvia three main channels: (1) Compton downscatter-ing of photons on cold electrons, (2) Coulomb colli- sions whereby nonthermal pairs lose their energy to coldplasma, and (3) photoionization.Below a few tens of keV, the dominant opacity sourceis photoionization. At these energies the photoelectronsare strongly coupled to thermal electrons; furthermore,the photoionization cross-section generally increases asthe radiation energy is reprocessed towards lower fre-quencies by repeated photoionizations and recombina-tions. For these reasons it is reasonable to assume thatall energy lost to photoionization is efficiently thermal-ized, i.e. occurs faster than any other timescale of inter-est.Between tens of keV and ∼
100 MeV most of thephoton energy loss is due to Compton downscattering.Above m e c ∼ ∼ MeV electron subsequently losesits energy either via IC scattering or Coulomb interac-tions with thermal electrons, both which typically occurfaster than the ejecta expansion time over which thepairs would experience adiabatic losses (Fig. 13). Inthe former case the photon is upscattered from the op-tical/UV to ∼ keV energies where it is thermalized byphoionization, while in the latter case the ∼ MeV elec-tron directly shares its energy with the thermal electronpool.Photon-matter pair production begins to dominatethe photon energy loss above a few tens of MeV. How-ever, only a fraction of the energy used for pair creationis eventually thermalized. The created secondary lep-tons have, on average, half of the energy of the primaryphoton. They cool mostly by IC and free-free emission,as well as by Coulomb collisions which become increas-ingly dominant as the lepton energy decreases (eqs. B11,B15). The thermalization efficiency of the lepton’s en-ergy depends on (1) the fraction of its energy lost toCoulomb collisions; (2) the shape and characteristic fre-quency of the pair’s IC and free-free radiation. Theenergy emitted in ∼ <
10 keV photons is photoelectricallyabsorbed and thermalized, while higher energy radia-tion either escapes, is lost by direct Compton (and ther-malized), and/or generates a subsequent generation ofelectron-positron pairs.The issue of thermalization is directly addressed byour Monte Carlo simulations in Section 3, which fol- Roughly given by the ratio of the lepton Lorentz factor belowwhich Coulomb losses dominate, γ (cid:63)(cid:63) (eq. B14), and the Lorentzfactor at which the lepton was created. Vurm & Metzger E (MeV) e ff E (MeV) e ff t (days)10 e ff t (days)10 e ff Figure 4.
Top Panels:
Effective optical depth through the supernova ejecta, τ eff , as a function of photon energy shown atdifferent snapshots in time after the explosion as marked. The calculations are performed for magnetar and ejecta parametersfit to the optical light curves of SN2015bn ( Left ) and SN2017egm (
Right ). We explore the dependence of the results on twoapproximations for the optical radiation field at late times, which serve as targets for γγ pair production. The cases shown assolid lines assume blackbody emission at the equilibrium temperature T eff = ( L opt / πσR ) / , where L opt is the optical/UVluminosity. The cases shown as dashed lines were instead calculated assuming a floor on the blackbody temperature at T eff = 4000K, as a crude proxy for the supernova nebular spectrum comprised of optical/NIR emission lines. We neglect nonthermal X-raytarget radiation in calculating the γγ optical depth, as is generally a good approximation at late times when τ γγ ∼ < BottomPanel:
Effective optical depth, now as a function of time after explosion for different photon energies as marked. low the diffusion of photons out of the spherical ejectacloud and self-consistently follow the electron-positronpair cascade initiated by the high-energy radiation. Tobriefly preface our results here (see Fig. 12 for details),for magnetized nebulae ( ε B >
0) Compton and photo-electric thermalization are found to contribute roughlyequally to the ejecta heating for the first few months af-ter explosion while the Thomson optical depth is high.However, photoelectric absorption comes to dominatethe ejecta heating at late times, once synchrotron emis- sion contributes a greater fraction of the nebula cool-ing. Coulomb heating also contributes appreciably dur-ing the transition period between the early Compton-and late photoelectric-dominated epochs.2.5.
Pair-loading of the pulsar wind and regulation ofthe mean energy of pairs entering the nebula
Consider a simple model in which the pulsar luminos-ity is dissipated in a localized radial zone, near or justprior to the termination shock (where the radial mo- amma-ray thermalization and leakage from millisecond magnetar nebulae Table 2.
Regimes of γγ pair creation regulation of the pulsar wind, in approximate temporal order of dominance (see AppendixD for details) High-energy γ -ray Target photon Regulated pair energy, γ † End time of regulation † IC IC ∼ / (2 √ θ ) ∼ − t ± , IC ≈ . ∼ ( m e c / θhν B ) / ∼ (10 − ) t / t ± , IC − syn ≈ ∼ /θ ∼ − t ± , th = min[ t γγ , t B ] ∼ . − ∼ ( m e c /hν B ) / ∼ (10 − ) t yr t ± , syn ∼ −
10 yr (eq. D45) † Numerical estimates assume fiducial parameters: B d ∼ G, v ej ∼ − , T eff ∼ − K ( θ ∼ − − − ). mentum flux of the wind matches the nebula pressure)and transferred to the pairs entering the nebula from theupstream wind region (Fig. 1; Kennel & Coroniti 1984;Sironi & Spitkovsky 2011). The mean energy gain perparticle, ∆ γ = L e ˙ N ± m e c , (33)depends on the ratio of the wind luminosity, L e , andthe number of pairs carried into the dissipation radius,˙ N ± . The latter is the sum of pairs injected directly fromthe engine on small scales (eq. 13) and those generated in situ via γγ pair creation (eq. 18). The high-energygamma-rays which create pairs in the wind can originateboth directly from the nebula (radiation from freshlyheated pairs) as well as from the supernova ejecta (in-wards scattering/diffusion following incomplete repro-cessing of gamma-rays from the nebula). This creates afeedback loop between ∆ γ and the number of generatedpairs in the wind: as the energy gain per particle in-creases, the resultant higher photon energies emitted bythe nebula increase pair production in the wind, whichin turn reduces the per-particle energy once those pairsreach the dissipation region.In a (quasi-)steady state the number M e of secondarypairs produced in the wind region per single leptonheated at the termination shock must be unity, i.e. M e ≈ . (34)If this was not the case (say, each heated lepton wereto instead generate more the one secondary), then ˙ N ± would increase and ∆ γ would correspondingly decrease.Typically, both the pair-production opacity and theoverall “gain” factor of the cycle decrease with decreas-ing ∆ γ , hence the system would readjust itself to attaina gain factor ∼ t ± , after which each pair-loading channel “breaks” and ∆ γ increases to the nextrelevant channel.At early stages, around and just after optical peak,the value of ˙ N ± due to γγ pair loading of the wind isso high that leptons entering the nebula do not generatesufficiently energetic photons to interact with the opti-cal/UV radiation field; the generation of secondary pairsat this stage is dominated by interactions of IC gammarays with other nonthermal photons (typically, X-raysand soft gamma rays also of IC origin). This phase gen-erally includes the epochs at which the ejecta becomesoptically thin to photon-matter pair production (Fig. 3),enabling gamma-ray emission in the Fermi
LAT rangeto escape. For sufficiently high values of the nebula mag-netization, the interaction of the IC-generated photonswith synchrotron targets can also dominate the windpair loading for times up to about a year.As the nebula expands and the compactness decreases,the value of ∆ γ increases until the cooling pairs radiatephotons of sufficient energy to also pair produce on the ∼ −
10 eV thermal optical photons. We explore thiscase here in some detail here, as it typically dominatesthe wind pair-loading on timescales of years after theexplosion, when the system is becoming transparent atthe highest gamma-ray energies. A more detailed ex-ploration of the results described below is provided inAppendix D.Figure 5 shows the numerically calculated multiplic-ity M e due to γγ pair creation in the pulsar wind bya single pair injected to the nebula. We use exactIC/synchrotron cooling rates and spectra assuming atarget radiation field dominated by thermal radiationof dimensionless temperature θ ≡ kT eff /m e c and showthe results as a function of the initial injected Lorentzfactor γ of the pair. Different lines show the results fordifferent values of the γγ optical depth τ γγ, th throughthe wind region and for different ratios of the nebula’smagnetic to thermal compactness, (cid:96) B /(cid:96) th (eq. 15).Unsurprisingly, Fig. 5 shows that the multiplicity in-creases with the initial particle energy γ , the mean en-ergy of the target radiation θ , and the gamma-ray op-2 Vurm & Metzger tical τ γγ, th . On the other hand, M e decreases with in-creasing (cid:96) B /(cid:96) th . This is because synchrotron photons(which dominate IC in cooling the pair for (cid:96) B /(cid:96) th ∼ > . −
1) have lower energies which are insufficient toproduce pairs on the optical photons (see eq. 26 andsurrounding discussion) .At early epochs t ± , IC ∼ < t ∼ < t B (eq. 17) when thenebula is still weakly magnetized ( (cid:96) B /(cid:96) th ∼ < . M e ∼ γ θ ∼ f ew . Themean Lorentz factor of particles entering the nebula,as well as the characteristic energy of the photons theyradiate (but not necessarily those that escape), willthus remain close to ∆ γ = γ ∼ /θ as long as thewind is opaque to gamma-rays above this energy (i.e. τ γγ, th > instead of theentire ejecta we see that the condition τ γγ, th > θ = kT eff / ( m e c ) ∼ − for T eff ≈ γ m e c ∼ m e c /θ ∼ γγ pair production at theseenergies ( t γγ ∼ yrs; eq. 12).On the other hand, Figure 5 also reveals that for highmagnetization (cid:96) B /(cid:96) th ∼ > . t ∼ > t B ; eq. 17) the self-regulation condition M e ∼ γ . It would then appear that the pairregulation process should break down at t ∼ t B evenif τ γγ, th > ε B is sufficiently large for t B (cid:28) t γγ thensynchrotron radiation dominates over thermal photonsas the targets for pair creation anyways (Table 2).Finally, at late times t (cid:29) t γγ once the wind is trans-parent to IC radiation, a final stage of regulated pairproduction can occur, this time by synchrotron gamma-ray photons on synchrotron X-ray targets. This regu-lation phase can last for a time t ± , syn ∼
10 yr, afterwhich the mean pair energy ∆ γ increases to the com-pletely “unregulated” value set by the intrinsic wind pairloading ( γ in ; eq. 18) or by the synchrotron burnoff limit( γ rad ; eq. 19). MONTE CARLO SIMULATIONS3.1.
Details of the Numerical Model
We perform three-dimensional Monte Carlo sim-ulations which track the energy evolution of elec-tron/positron pairs and photons in the magnetar nebula The thermal compactness of the ejecta (cid:96) th must be replaced bythe somewhat higher compactness of the wind region, (cid:96) th , n ≈ ( R ej /R n ) (cid:96) th arising due to its smaller radius R n ≈ R ej / e l B / l th = 0 l B / l th = 10 l B / l th = 10 = 10 = 10 = 10 ,th = 1 ,th = 0.3 Figure 5.
Pair multiplicity M e generated in the pulsarwind by γγ interactions as a function of 4 γ θ , where γ isthe average Lorentz factor of heated leptons entering thenebula and θ ≡ kT eff /m e c is the dimensionless temperatureof the background target radiation field (assumed here to bethermal radiation with T eff ∼ − K, correspondingto θ ∼ − − − ). Separate lines show the results fordifferent values of τ γγ, th (the γγ pair creation optical depththrough the wind region; eq. 30) and the ratio of magneticto thermal compactness (cid:96) B /(cid:96) th in the nebula (eq. 16). Asdescribed in the text, γγ pair loading of the pulsar windregulates M e ∼ γ ∼ /θ so long as τ γγ, th ∼ > (cid:96) B /(cid:96) th ∼ < .
1, conditions typically satisfied for months toyears following the supernova explosion.
Table 3.
Magnetar model parameters used in Monte CarlosimulationsSN B d P M ej v ej κ opt t e - (10 G) (ms) ( M (cid:12) ) (km s − ) (cm g − ) (d)2015bn 0.43 2.6 11.7 5400 0.1 542017egm † † Broadly motivated by Bayesian fits of the Kasen &Bildsten (2010) model to supernova optical light curve data(Nicholl et al. 2017a). and supernova ejecta, as well as the outwards radial dif-fusion of both high-energy and optical photons throughthe ejecta. We employ a two-zone approach, which con-siders separately a homogeneous inner “wind/nebula”zone and an outer “ejecta” zone (Fig. 1). Photons andelectron/positron pairs can occupy both zones, but ionsand their associated electrons are restricted to the ejectashell. Conversely, the magnetic field is restricted to the amma-ray thermalization and leakage from millisecond magnetar nebulae R n = R ej / R ej and zero outside thisrange. The time-dependent Thomson optical depth ofthe ejecta shell τ T ∝ t − (eq. 2) is determined from theassumed fixed ejecta velocity v ej and mass M ej . We donot account for acceleration of the ejecta by the nebula,because throughout most of the evolution of interest thenebula/ejecta is radiative (such that most of the engine’sluminosity is lost to radiation rather than converted intobulk kinetic energy). The true ejecta distribution willundoubtedly be more complex than the constant densityshell we have assumed, for instance being compressedinto a thin radial shell by the nebula (e.g. Kasen & Bild-sten 2010) or characterized by inhomogeneities and mix-ing resulting from Rayleigh-Taylor instabilities at thenebula-ejecta interface (e.g. Chen et al. 2016; Suzuki &Maeda 2020). However, we do not expect our qualitativeconclusions to depend on such details.Primary pairs are injected into the nebula at the ra-dius R inj = R n corresponding to the termination shockat the rate ˙ N ± specified by Equation (13) for an assumedGoldreich-Julian pair multiplicity µ ± = 10 . The mul-tiplicity is uncertain theoretically, but our key resultsturn out to be insensitive to its precise value becausethe mass loading of the wind over timescales relevantto optical and gamma-ray observations is dominated bysecondary pairs generated by γγ interactions in the windzone (Section 2.5). The latter are explicitly tracked af-ter being generated as they are advected outwards in thewind.Upon reaching the termination shock radius, both theprimary and secondary pairs are energized based on thecurrent pulsar spindown luminosity, L e (eq. 4). Depend-ing on the total rate that particles enter the shock, ˙ N ± ,each particle experiences an energy gain of ∆ γ ∝ L e / ˙ N ± (eq. 33). All pairs entering the shock receive the sameenergy gain, regardless of their initial energy and we donot account for potential nonthermal particle accelera-tion at the shock. At early times, when γγ pair-loadingdominates the particle budget of the wind, the two-way feedback between the number of secondary pairsgenerated in the wind zone and the dissipated energyper particle ∆ γ (Section 2.5), introduces oscillations in∆ γ causing some “jitter” in our results on top of thatpresent due to the Monte Carlo sampling.We assume that pairs radiate their energy in the neb-ula isotropically with negligible bulk velocity. In prin-ciple this assumption might be questioned because atearly times of interest the γγ pair creation − and con-comitant mass-loading and heating of the wind − occur over a region of finite radial thickness ahead of the windtermination shock roughly equal to the mean-free pathof high-energy photons to γγ interactions. If the hotpairs entrained in the wind were to cool radiatively overthis radial scale faster than being advected into the neb-ula, then the engine’s luminosity would instead be radi-ally beamed due to the relativistic motion of the wind,violating the assumption of isotropic emission and alter-ing the energy spectrum of the pairs’ emission.In Appendix C we show that IC cooling can be ne-glected within the pair loading region so long as thepulsar wind remains relativistic all the way to the termi-nation shock. The condition for negligible synchrotroncooling is somewhat more stringent and requires eitherthe wind possess a bulk Γ (cid:29) f − / ∼ −
100 or lowmagnetization σ (cid:28) f th Γ by the radius of pair-loading.These conditions could be satisfied if the Poynting fluxof the wind is largely converted into heat or bulk kineticenergy by the termination shock (e.g. as is inferred tooccur in the Crab pulsar wind; Kennel & Coroniti 1984),or if the process of γγ pair-loading the wind itself actsto reduce its magnetization. Though beyond the scopeof this paper, we plan to explore the latter possibility infuture work (see Section 4.1 for additional discussion).As long as the pairs do not cool in the relativistic windzone, the precise location at which they are energized isnot critical because the particles radiate most of theirenergy isotropically once reaching the nebula. This al-lows us to forgo treating the wind zone as a separateregion regarding synchrotron and IC processes.The magnetic field within the wind/nebula region isgiven by Equation (14), i.e. we assume the magneticenergy to be a fixed fraction ε B of the total energy de-posited by the engine over the ejecta expansion time. Wegenerally assume modest values of the nebula magneti-zation ε B ∼ <
1, motivated by the likelihood of magneticreconnection within the sub-relativistically expandingnebula (e.g. Begelman 1998; Porth et al. 2013) and thepotential reduction in wind magnetization even prior tothe termination shock due to γγ pair-loading as men-tioned above.All of the photon destruction and creation processesdiscussed in Section 2.3 are included in both the neb-ula and ejecta. In addition to the high-energy (non-thermal) radiation, we also explicitly follow the creationand transport of optical (thermal) photons. The opti-cal photons are created by thermalization processes de-scribed in Section 2.4, namely Compton downscatteringon cold electrons, Coulomb losses by nonthermal elec-trons/positrons and photoionization. While the detailedspectrum of optical radiation is not considered, the codefollows the frequency-integrated energy of the optical ra-4 Vurm & Metzger diation field represented by MC particles. Their inter-actions with the ejecta material as they diffuse outwardare assumed to take place with a constant opacity κ opt .The radial expansion of the ejecta is taken into accountin the interaction events (modeled as coherent scatteringin the local instantaneous rest frame), thereby automat-ically accounting for adiabatic losses and hence suppres-sion of the early optical radiation when τ opt > c/v ej .As summarized in Table 3, we adopt parameters forthe ejecta ( M ej ; v ej ; κ opt ) and magnetar engine ( B d ; P )which are broadly motivated by model fits to the opticallight curve data on two well-studied SLSNe, 2015bn and2017egm (Nicholl et al. 2017a, 2018). Once the basicejecta and engine parameters have been specified, themain free parameter is the nebula magnetization, whosevalue we vary from ε B = 0 to ε B = 10 − (similar to theCrab Nebula). 3.2. Results
Figures 6 – 11 summarize our results for for SN2015bnand 2017egm for different values of the nebula magne-tization. For the illustration of key observables, suchas the optical and gamma-ray light curves, we presentresults for both supernovae. However, given that thetwo models exhibit a qualitatively similar evolution, forthe more interpretation-focused figures we present justthe SN2015bn cases. Insofar as these two supernovaeare also fairly representative of the larger SLSNe popu-lation (Fig. 2), we expect qualitatively similar results tothose presented to hold more broadly.The optical light curves (Figs. 6, 7) around and justafter peak are largely insensitive to magnetization and,in fact, more generally to the details of energy depo-sition and reprocessing of high-energy radiation withinthe nebula/ejecta. This is because reprocessing at thisearly stage is nearly 100% efficient; thus, most of thehigh-energy radiation gets downgraded into the opticalband where opacity is lower compared to higher frequen-cies. Defining a “thermalization” optical depth, τ therm = − ln(1 − L opt /L e ) , (35)we find τ therm > t ∼ <
100 days (Figure 10).The opacity κ th ≡ ( τ therm /τ T ) κ es which corresponds to τ thermal is equivalent to the gamma-ray thermalizationopacity introduced in previous SLSNe modeling (e.g.,Wang et al. 2015; Chen et al. 2015), with the impor-tant difference that its value is calculated here self-consistently and is not constant in time.Our theoretically predicted optical light curves risesomewhat more abruptly than SLSNe observations(shown as black points in Figs. 6, 7) prior to opticalmaximum. This is partly a consequence of our assump-tion of initially cold ejecta and the sole energy source being the centrally concentrated nebula. In the physicalcase, other sources of volumetrically-distributed heat-ing of the ejecta not accounted for in our model maycontribute to the thermal emission at early times andflatten the rise. Such additional heating sources includethe shock driven into the ejecta by the nebula at earlytimes (e.g. Kasen et al. 2016; Chen et al. 2016; Suzuki &Maeda 2019) or a relativistic jet from the central engine(e.g., Margalit et al. 2018b).Starting a few months after the explosion, τ therm fallsbelow unity and high-energy radiation begins leakingout of the ejecta. Most of the gamma-rays which es-cape at early times are generated by IC scattering in thenebula. The ejecta first become transparent in the en-ergy range 0 . −
10 GeV, within the bandpass accessibleto
Fermi -LAT (Section 2.3). The corresponding ∼ GeVlight curves (red lines in Figs. 6, 7) rise to a maximumat ∼
200 days and t ∼
100 days, respectively, at whichstage they carry a sizable fraction of the engine lumi-nosity; any higher energy gamma-ray emission at thisepoch is also reprocessed down into the ∼ GeV band.The emission properties at yet higher energies ∼ >
100 GeV rely on both the pulsar wind’s ability to in-ject sufficiently energetic pairs to create such photons,as well as the ability of those photons to escape the neb-ula and ejecta in face of γγ absorption. Regarding thefirst condition, Figure 11 shows that the average energyof the injected pairs has already increased to ∆ γ ∼ > ( ∼ >
100 GeV) by ∼
200 days (as set by the regulated γγ pair-loading in the wind; Section 2.5). These energeticpairs mainly cool by IC radiation on optical targets inthe Klein-Nishina regime (at the higher magnetizations,synchrotron also competes), depositing most of their en-ergy into photons with energies comparable to the cool-ing lepton (see eq. 24).The delayed onset of ∼ >
100 GeV until t ∼ −
200 days is thus mainly controlled by opacity, whichat these energies is dominated by γγ interactions on theoptical radiation field (Fig. 4). Because opacity controlsthe ability of the highest-energy gamma-rays to escape,the peak timescale and luminosity of the ∼ >
100 GeVemission is strongly dependent on the thermalization ef-ficiency (eq. 35) insofar as the latter determines the frac-tion the engine’s luminosity which is reprocessed into op-tical target radiation. A steeper(shallower) decay of theoptical light curve hastens(lengthens) the γγ optically-thin transition and generally increases(decreases) thepeak luminosity in the ∼ >
100 GeV bands.On the other hand, a very steeply decaying opticallight curve (e.g. the ε B = 0 cases in Figs. 6, 7) alsoacts to suppress the late high energy emission. This isbecause, absent synchrotron cooling, a weak target ra- amma-ray thermalization and leakage from millisecond magnetar nebulae t (days) L ( e r g s ) B = 0 Optical0.1 - 500 GeV0.2 - 30 TeV L e t (days) L ( e r g s ) B = 10 Optical0.1 - 500 GeV0.2 - 30 TeV L e t (days) L ( e r g s ) B = 10 Optical0.1 - 500 GeV0.2 - 30 TeV L e t (days) L ( e r g s ) B = 10 Optical0.1 - 500 GeV0.2 - 30 TeV L e Figure 6.
Optical and gamma-ray light curves of magnetar-powered supernovae, calculated for different values of the nebulamagnetization ε B = 0 , − , − , − as marked above each plot. The magnetar and ejecta parameters (Table 3) follow modelfits to the optical light curve data for SN2015bn (shown as black dots; Nicholl et al. 2016a). A dot-dashed line shows the engineluminosity L e (eq. 6). diation field for IC cooling hastens the transition of thenebula into an adiabatic regime, after which the injectedpairs are no longer able to cool efficiently over a dynam-ical time. This results in a turnover in the light curvesat all energy bands for ε B ≈ ε B = 10 − and ε B = 10 − models. Gamma-ray emission initiallyappears in the 0 . −
10 MeV range around peak opticallight. Before t ≈
100 days the emission is significantlyattenuated by Compton recoil losses. At the low-energyend the spectrum is truncated by bound-free absorption,whereas the high-energy turnover represents the maxi-mum energy of IC upscattered optical photons. The lat-ter is limited by the fact that the energies of the pairs in the nebula responsible for the upscattering at this stageare regulated to a low value (Figure 11) due to efficientpair production on nonthermal radiation (Table 2).As the transparency window expands, the spectrumbroadens and becomes limited by photo-electric absorp-tion and thermal pair production at the low and high-energy ends, respectively. At t ∼
300 days the spectrumis relatively flat in νF ν space, softened relative to theclassical cooling spectrum νF ν ∝ ν / by pair-photoncascades. As the overall compactness decreases, the cas-cade ceases and the spectrum correspondingly hardens.The bulk of the gamma-ray energy at ∼ Vurm & Metzger t (days) L ( e r g s ) B = 0 Optical0.1 - 500 GeV0.2 - 30 TeV L e t (days) L ( e r g s ) B = 10 Optical0.1 - 500 GeV0.2 - 30 TeV L e t (days) L ( e r g s ) B = 10 Optical0.1 - 500 GeV0.2 - 30 TeV L e t (days) L ( e r g s ) B = 10 Optical0.1 - 500 GeV0.2 - 30 TeV L e Figure 7.
Same as Figure 7 but for magnetar model parameters fit to SN2017egm (optical light curve data from Nicholl et al.2017a; Bose et al. 2018).
Starting approximately 2 − γγ optical depth drops below unity and all the ra-diated TeV photons can escape from the source. Thisalso marks the time when the γγ pair-loading of thewind begins to wane, causing the energy per lepton toincrease (Fig. 11; see discussion in Section 2.5). As aresult, the IC scattering evolves even deeper into theKlein-Nishina regime, which has two significant conse-quences: 1. synchrotron emission can become the dom-inant cooling mechanism even if (cid:96) B /(cid:96) th (cid:28)
1, and 2. incase the electron cooling is still dominated by IC radia-tion (e.g. in the ε B = 10 − model) then the transitionto the adiabatic regime takes place earlier than it wouldfor lower ∆ γ ; both have a suppressing effect on high-energy gamma ray emission. The emergence of an en-ergetically dominant synchrotron component is seen at t ∼
600 and 2000 days in the ε B = 10 − and 10 − cases,respectively, at 10 −
100 keV (magenta and green lines,right panels of Fig. 8). Prior to this time most of theescaping gamma-rays were generated by IC scattering.The time-integrated energy radiated in the synchrotroncomponent is seen to be greater for larger values of ε B .Figure 10 shows the time evolution of the thermal-ization optical depth of the ejecta (eq. 35), shown sep-arately on its own (left panel) and normalized to theThomson optical depth τ T (right panel). The τ therm /τ T ratio is proportional to the thermalization opacity κ th employed in previous SLSNe model and shown on theright vertical axis. The effective κ th is larger for highervalues of ε B due to the greater role of synchrotron emis-sion cooling in the nebula than IC scattering for strongermagnetic fields: the lower energies of the synchrotron amma-ray thermalization and leakage from millisecond magnetar nebulae h (GeV) E ( e r g ) Cumulative spectra ( B = 10 ) t = 100 days t = 200 days t = 300 days t = 600 days t = 1000 days t = 2000 days h (GeV) L ( e r g s ) NuSTAR Fermi -LAT
Veritas / MAGIC / CTA
Instantaneous spectra ( B = 10 )10 h (GeV) E ( e r g ) Cumulative spectra ( B = 10 ) t = 100 days t = 200 days t = 300 days t = 600 days t = 1000 days t = 2000 days h (GeV) L ( e r g s ) Instantaneous spectra ( B = 10 ) t = 100 days t = 200 days t = 300 days t = 600 days t = 1000 days t = 2000 days Figure 8.
Cumulative and instantaneous energy spectra of the escaping high-energy radiation from the ε B = 10 − and 10 − models for SN2015bn, corresponding to the second and last panels of Fig. 6. photons result in them being thermalized with higher ef-ficiency. The rise of τ therm /τ T at hundreds of days occursonce synchrotron emission begins to dominate the ther-malized luminosity (even if synchrotron does not nec-essarily yet dominate the total nebula cooling because (cid:96) th (cid:29) (cid:96) B ; Fig. 9). The abrupt plunge in τ therm at t ∼ ε B = 0) case is again theresult of the nebula becoming adiabatic to IC cooling asa result of the reduced seed field due to energy escapingthe nebula. All models depart at late times from the ap-proximately constant value of τ therm /τ T ≈ . γγ pair pro-duction in the pulsar wind on IC-generated X-rays reg-ulates the injected pair energy to a relatively low value(∆ γ ∼ / √ θ ∼ Vurm & Metzger t (days) , T injth,n ; B = 0 th,n ; B = 10 ; B = 10 ; B = 10 ; B = 10 ; B = 10 ; B = 10 Figure 9.
Time evolution of injection, thermal andmagnetic compactnesses and Thomson optical depth(eqs. 8,10,15,2) for the SN2015bn models shown in Figure6.
The channels for loading secondary pairs into the pul-sar wind on nonthermal IC radiation becomes less ef-ficient a few hundred days after the explosion. Theavailable energy per lepton then increases until thermalpair production takes control of the pair loading, en-forcing ∆ γ ∼ /θ (Table 2). The overall thermalizationefficiency decreases, as the bulk of the dissipated lumi-nosity now escapes as MeV-GeV photons. Albeit withdecreasing efficiency, the low-energy part of the IC spec-trum continues to contribute to thermalization via thephotoelectric and Compton recoil processes (Figure 12),roughly according to the effective opacity shown in Fig-ure 4. In addition, the lower-energy ( γ (cid:28) ∆ γ ) electron-positron pairs generated within the ejecta by photon-matter and photon-photon pair production contributeto heating/thermalization via Coulomb losses.However, after ∼ −
200 days synchrotron emissionbecomes the dominant contributor to the thermalized lu-minosity, unless the nebula magnetization is extremelylow. At this stage, the synchrotron spectrum typicallypeaks in the X-ray domain (eq. 26) and a large frac-tion of its luminosity is thermalized by photo-electricabsorption. In this regime one can derive a simple ap-proximate relation between the injected and thermalizedluminosities (i.e. thermalization efficiency). Assumingthat the nebula pairs cool predominantly by IC upscat-tering thermal photons and the cooling is still fast com-pared with the dynamical time ( t < t c ; eq. 23), one can write L opt ≈ L e ˙ γ syn ˙ γ IC (1 − e − τ bf ) = L e u B , n F KN u th , n (1 − e − τ bf ) , (36)where u B , n and u th , n are magnetic and thermal radi-ation energy densities within the nebula, respectively,and τ bf is an appropriately defined “average” photo-electric opacity of the synchrotron photons. Using equa-tion (14), the magnetic energy density can be expressedas u B = B π = L e πcR ε B cv n , (37)whereas the thermal energy density is u th = L opt / (4 πcR ). Substituting u B and u th into equation(36), we obtain τ therm ≈ L opt L e ≈ (cid:114) ε B cF KN v n (1 − e − τ bf ) . (38)For example, in our SN2015bn model with ε B =10 − , synchrotron emission is efficiently absorbed the bybound-free process for the first ∼ τ bf (cid:29) t ≈ γ ∼ /θ , such that F KN ≈ .
05 (Sec-tion 2.5). This, along with v n ≈ v ej / ≈ − ,yields τ therm ≈ L opt /L e ≈ . τ therm (Figure 10, left panel), owing to contributionsfrom Compton and Coulomb thermalization at earliertimes as well as somewhat higher ∆ γ (and thus lower F KN ) towards the end of the considered time interval.These are also the reasons for the slight upward curva-ture of τ therm , which would otherwise remain constantas long as the above assumptions are valid. DISCUSSION4.1.
Late-time optical light curves as probes of thenebula magnetization
If SLSNe are powered from within by reprocessing ofa central gamma-ray source, this should manifest in-directly in the late-time behavior of their optical lightcurves. Our self-consistent models of the nebula gamma-ray emission and its thermalization provide an opportu-nity to explore this question on solid theoretical footingand explore its implications for the magnetar model.Nicholl et al. (2018) fit SLSNe light curves to themagnetar model, including a suppression factor f th =(1 − e − τ therm ) on the engine luminosity to account forincomplete thermalization of the engine gamma-rays atlate times, where τ therm ∝ κ th τ T ∝ κ th t − and κ th is amma-ray thermalization and leakage from millisecond magnetar nebulae t (days) t h e r m ; T therm : B = 0 therm : B = 10 : B = 10 : B = 10 t (days) t h e r m / T B = 0 B = 10 B = 10 B = 10 radioactive decay t h ( c m g ) Figure 10.
Evolution of the thermalization optical depth τ therm (eq. 35) for the SN2015bn models in Figure 6. On the rightpanel we have normalized τ therm to the Thomson optical depth τ T ∝ t − (eq. 2), which is shown separately with a dashed lineon the left panel. The ratio τ therm /τ T = κ th /κ es can also be interpreted in terms of the gamma-ray thermalization opacity κ th introduced in previous SLSNe modeling (e.g. Wang et al. 2015; Nicholl et al. 2017c; which however assume a constant value of κ th ). A dotted horizontal line shows the constant κ th = 0 . κ es which approximates thermalization of Ni and Co radioactivedecay in ordinary hydrogen-poor supernovae (Swartz et al. 1995). t (days) B = 0 B = 10 B = 10 B = 10 Figure 11.
Energy received per lepton at the wind termi-nation shock for the SN2015bn models (Fig. 6). Shown forcomparison are: γ in (eq. 18), the mean energy pairs wouldacquire entering the nebula absent γγ loading; 1 /θ , an ana-lytic estimate of the self-regulated particle energy when γγ pair production on thermal radiation regulates the wind pairloading (see Section 2.5, Fig. 5); and γ rad (eq. 19), the maxi-mum particle energy corresponding to the synchrotron burn-off limit. taken to be a constant gamma-ray thermalization opac-ity (Wang et al. 2015; Chen et al. 2015). At late times,when τ therm (cid:28)
1, this model predicts a decay of the su-pernova luminosity L opt ∝ L e f th ∝ L e τ therm ∝ t − for amagnetar engine with L e ∝ t − (eq. 4), in reasonable ac- cord with the best-fit power-law decay L opt ∝ t − . mea-sured for SN2015bn over the time interval t = 200 − κ th ∼ . − . g − inferredfrom the SLSNe light curve fits (e.g. Nicholl et al. 2017c)are similar to the value κ th ∼ .
02 cm g − expectedfrom the thermalization of ∼ MeV gamma-rays gener-ated by the radioactive decay of Ni → Co → Fe inordinary stripped envelope supernovae. In the case ofradioactive decay, κ th is indeed expected to be roughlyconstant over the first ∼ ∼ MeV spectrum of ra-dioactive decay products (e.g., Fig. 8) and there is noreason a priori to expect a similar value of κ th to or-dinary supernovae, much less one which is constant intime.On timescales of months after the explosion our MonteCarlo simulations do in fact reveal an effective value τ therm /τ T ∼ . − κ th ∼ . − . g − (Fig. 10). However, the predicted time evolu-tion strongly deviates from κ th = const , with κ th ini-tially decreasing before reaching a minimum and thenbeginning to rise. At late times the reprocessed opticalluminosity approaches a fixed fraction of the engine lu-minosity ( τ therm /τ T ∝ κ th ∝ t ; eq. 38). The transition0 Vurm & Metzger t (days) L ( e r g s ) B = 0 L th L th : Compton L th : Photoelectric L th : Coulomb L opt L e t (days) L ( e r g s ) B = 10 L th L th : Compton L th : Photoelectric L th : Coulomb L opt L e t (days) L ( e r g s ) B = 10 L th L th : Compton L th : Photoelectric L th : Coulomb L opt L e t (days) L ( e r g s ) B = 10 L th L th : Compton L th : Photoelectric L th : Coulomb L opt L e Figure 12.
Thermal luminosity L th responsible for powering the optical light curves of SLSNe, and its partition into differentchannels of thermalization (Compton, Coulomb, Photoelectric) described in Section 2.4, for the models of SN2015bn shown inFigure 6. Photoelectric absorption of soft synchrotron X-ray photons dominate the thermalization of the nebula radiation atlate times, unless the nebula magnetization is very low (see eq. 38 and surrounding discussion). to a rising κ th is driven by synchrotron emission be-coming the dominant mode of thermalization and henceoccurs earlier for larger values of the nebula magnetiza-tion, ε B . In particular, to match the steep optical lightcurve decay L opt ∝ t − . of SN2015bn out to t ∼ > ε B ∼ < − (Fig. 6). For SN2017egm,matching the optical luminosity at t ∼
250 d requires ε B ∼ < − (Fig. 7).Stated another way, our models with ε B ∼ > − − − significantly overpredict the late-time optical lightcurves of SN2015bn and SN2017egm. If the magnetarscenario for these SLSNe (as commonly understood) iscorrect, then our results show that highly efficient dis-sipation of the wind’s magnetic field must occur in the wind or nebula (near the light cylinder, the energy of thepulsar wind is carried almost entirely in Poynting flux;Goldreich & Julian 1970). Magnetic dissipation can takeplace in the nebula due to non-axisymmetric instabili-ties arising from the dominant toroidal magnetic fieldgeometry (Begelman 1998); however, the relatively highresidual magnetization observed in ordinary pulsar windnebulae such as the Crab Nebula ε B ∼ > − (e.g. Ken-nel & Coroniti 1984; Begelman & Li 1992) suggest thatsuch a dissipation process alone may not be sufficient toreach the low values of ε B required to explain SLSNe.One mechanism for converting Poynting flux into ki-netic or thermal energy is magnetic reconnection in thepulsar wind ahead of the termination shock, as may oc- amma-ray thermalization and leakage from millisecond magnetar nebulae γγ interactions (Section 2.5).Loading of the wind with relativistically hot pairs maynot only act to decelerate the outflow, but also reducethe wind magnetization as a result of the electric currentinduced by the deposited pairs prior to their isotropiza-tion in the co-moving frame of the wind (D. Giannios,private communication). Though beyond the scope ofthe present study, we plan to explore this possibility infuture work.Absent a viable mechanism for dissipating the neb-ula’s magnetic field with extraordinary efficiency, weconclude that the magnetar model for SLSNe aspresently envisioned may be theoretically challenged.We encourage additional late-time optical observationsof SLSNe to ascertain whether the steepening seen inSN2015bn and SN2017egm are generic and to searchfor the predicted late-time flattening arising from theturnover in τ therm /τ T .4.2. Prospects for gamma-ray detection of SLSNe
Our models predict that high-energy gamma-ray emis-sion will accompany SLSNe months to years after theexplosion, starting at ∼ MeV energies ( t ∼ >
50 days) andthen moving up to ∼ GeV ( t ∼ >
100 days) and eventually ∼ TeV energies ( t ∼ > . −
600 GeV luminosi-ties of SLSNe by
Fermi
LAT were compiled by Renault-Tinacci et al. (2018). On an individual basis these limitson Type I SLSNe, L LAT ∼ < − erg s − over a6-month window following the explosion, are in generalnot constraining on the luminosities predicted by ourmodels of ∼ < few 10 erg s − (Figs. 6, 7, 8).Higher energy gamma-rays ∼ >
100 GeV can be ob-served by atmospheric Cherenkov telescopes, such as thefuture planned Cherenkov Telescope Array (CherenkovTelescope Array Consortium et al. 2019), which can per-form pointed observations that achieve deeper νF ν sen-sitivity than obtained by Fermi
LAT in survey mode. However, these advantages are mitigated by the fact thatthe ∼ TeV gamma-ray emission predicted by our modelsonly rises on a timescale of years, due primarily to thehigh γγ opacity of the supernova optical light (Fig. 2).Our low- ε B models, which best reproduce the opticallight curves of SLSNe, fortunately also possess the great-est TeV luminosities (Figs. 6,7). However, even in thesebest-case scenarios the predicted ∼ TeV peak luminosi-ties are an order of magnitude smaller than that in the ∼ GeV band and a factor of several below the instanta-neous engine luminosity at this late epoch.Despite these challenges, observational efforts to ob-serve nearby and optically bright SLSNe in the gamma-ray bands on timescales of months to years are stronglyencouraged. A gamma-ray detection would representa convincing confirmation of engine-powered models(Kotera et al. 2013; Murase et al. 2015).4.3.
Implications for late-time spectral features
Nicholl et al. (2016b) propose that the atypical oxygenemission line features observed in the nebular spectra ofSN2015bn may be formed at the dense inner edge ofthe ejecta shell and arise due to excitation by radiationor energetic particles from the central engine (see alsoJerkstrand et al. 2017; Nicholl et al. 2019). Milisavljevicet al. (2018) argue that similar features in the nebularspectra of the (atypical but non super-luminous) TypeIb SN 2012au result from photoionization of the ejectashell by a central pulsar wind nebula.Our models explicitly identify what physical processesare responsible for heating the ejecta in engine-poweredmodels. As shown in Figure 10, at early times t ∼ <
200 da sizable fraction of the optical light curve is poweredby Compton thermalization, in which the upscatteredelectron shares a portion of its energy directly with theplasma via Coulomb scattering and another part viaphotoionization from secondary photons created by ICscattering off the optical radiation. By contrast, at latetimes t ∼ >
200 d, almost all of the thermalization is dueto photoionization by synchrotron photons. Thus, inso-far as “cosmic ray” versus photo-ionization energy de-position will affect the ionization state evolution of theejecta in distinct ways, our models’ predictions for thetemporal and radial distribution of the heating throughthese different channels, could serve as important inputto future models of SLSNe nebular spectra.4.4.
Late-time synchrotron radio emission
As electrons cool in the nebula via synchrotron radi-ation, they will produce emission extending down intothe radio bands ∼ <
100 GHz (Murase et al. 2016; Met-zger et al. 2017; Omand et al. 2018). The supernova2
Vurm & Metzger ejecta typically become optically thin at GHz frequen-cies on a timescale of several decades after the ex-plosion (e.g. Margalit et al. 2018a). Eftekhari et al.(2019) discovered radio emission from the location ofSLSN PTF10hgi at t ≈ . γγ pair deposition in the wind has typically ended( t ∼ > t ± ). Therefore, the mean energy per pair injectedinto the nebula at these lates time can be very large∆ γ ∼ min[ γ in , γ rad ] ∼ − (eqs. 18,19) for phys-ical values of the primary Goldreich-Julian pair mul-tiplicity of the wind µ ± ∼ < (Timokhin & Harding2019). Given that these high average energies result insynchrotron emission peaking in the X-ray or gamma-ray band (eq. 26), we expect that pulsar/magnetar windnebulae may be challenged to produce the radio lumi-nosity seen in PTF10hgi, which is a significant fractionof the total magnetar spin-down power on this timescale.Other sources of electron injection and heating, suchas those which accompany mildly relativistic ejectionsof baryon-rich matter from magnetically powered flaresand which can shock heat electrons entering the nebulato Lorentz factors ∼
100 in the appropriate range forradio emission, are potentially more promising for pro-ducing luminous late-time radio emission from magnetarengines (Margalit & Metzger 2018).4.5.
Application to other engine powered transients
Although we have focused on engine and ejecta param-eters appropriate to SLSNe, the general physical set-uppresented may apply to other optical transients whichhave been considered to be powered by a millisecondpulsar/magnetar or accreting black hole central engine.The primary requirement to preserve the qualitative pic-ture is only that the engine generate a relativistic pairoutflow behind an expanding ejecta shell, with a suffi-ciently high per-particle energy to activate the γγ regu-lation process described in Section 2.5.As one prominent example, the optical and X-ray/ γ -ray light curves of FBOTs such as AT2018cow may bepowered by a magnetar or black-hole formed in a corecollapse explosion characterized by ejecta significantlyless massive ( M ej ∼ < M (cid:12) ) and faster ( v ej ∼ > . M ej ∼ < . M (cid:12) ; v ej ∼ > . c ) thanin the FBOT case. Yet another potential application ofour model is to tidal disruption events (TDE) of stars bymassive black holes (e.g., Rees 1988), in which gravita-tional or rotational energy released as the stellar debrisis accreted by the supermassive black hole may be “re-processed” by the unbound debris from the disruption,playing an important role in powering the optical lightcurves of TDEs (e.g., Guillochon et al. 2014; Metzger &Stone 2016)Throughout Section 2 we have attempted couch keyphysical quantities in terms of dimensionless parame-ters such as compactness, which can be readily scaledto other systems. The peak time of optical transients isusually given by the photon diffusion time (eq. 7), whichfor fixed optical opacity scales as t pk ∝ M / v − / . Fora magnetar engine with a spin-down time t e ∼ < t pk , thepeak optical luminosity from Arnett’s Law is L pk ∼ L e ( t pk ) ∝ B − t − ∝ B − M − v ej (eq. 4) and hence B d ∝ M − / v / L − / .With this information in hand, we can now scale theother key timescales in the problem (Table 1) to thepeak time t pk and luminosity L pk . This gives: t γγ t pk ∝ M − / v − / L / ; t T t pk ∝ v − / (39) t B t pk ∝ constant ; t c t pk ∝ L / M − / v − / (40)These ratios reveal that, even scaling to their fasteroptical rise times, FBOTs and neutron star mergertransients with lower M ej and higher v ej will generallypass through the critical stages of evolution faster thanSLSNe. In particular, their shorter time spent at highgamma-ray compactness (smaller t γγ ) could lead to theinjected energy of the nebular pairs ∆ γ (and their re-sulting gamma-ray emission) rising more rapidly in time.Given also that the ejecta becomes Thomson thin faster(smaller t T ), the engine’s luminosity may begin to leakout as gamma-rays faster, causing the reprocessed opti-cal emission after peak light to decrease faster. Further-more, our implicit assumption throughout this work has amma-ray thermalization and leakage from millisecond magnetar nebulae ∝ M ej /v .The X-ray/gamma-ray emission observed over the firstseveral weeks to months of AT2018cow could be at-tributed to synchrotron emission from a magnetar orblack hole nebula escaping from a highly-ionized aspher-ical ejecta shell (e.g. Margutti et al. 2019). The engine-powered X-ray light curve furthermore showed a breakfrom L X ∝ t − to ∝ t − at approximately 30 days af-ter the explosion, potentially consistent with the spec-tral energy distribution of the nebula shifting to higherphoton energies once γγ regulation of the engine’s windmass-loading ceases. In SLSNe this transition occursroughly on the timescale t γγ ∼ t pk ∼ ∼ Caveats and uncertainties
The model we have presented is subject to several un-certainties and simplifying assumptions which requireclarification in future work. Some of these assumptions,such as of isotropic emission by the energized pairs afterentering the nebula, were already discussed in Section3.1 (see also Appendix C).Another earlier-mentioned uncertainty is the intrinsicGoldreich-Julian pair multiplicity of young pulsar winds, µ ± . However, unless the multiplicity is orders of magni-tude larger than generally assumed in pulsar wind neb-ulae, our results are not dependent on µ ± because thewind mass-loading is dominated by secondary pairs gen-erated by γγ interactions in the wind just ahead of thetermination shock (Fig. 11). For the same reason, theconclusions we have drawn for relativistic pulsar windswould be similar if the engine were instead the ultra-relativistic jet of an accreting black hole carrying thesame power and a similarly high initial per-particle en-ergy (though not necessarily if the engine were a mildlyrelativistic accretion disk wind with a comparatively lowper-particle energy).Our calculations neglect any dynamical effects of cool-ing in the nebula. Except in the case of an unmagne-tized nebula, the pairs remain radiative for decades tocenturies after the explosion (eq. 23). Absent additionalpressure support (e.g. turbulence or an ion componentof the magnetar ejecta), the ram pressure of pulsar windcould compress the nebula into a thin layer, of radialthickness ∆ R n /R n ∼ v ej /c ∼ .
05, potentially leading to dynamical effects or instabilities not captured by ourmodel.More broadly, our calculations assume the nebula andejecta shell are spherically symmetric. Suzuki & Maeda(2020) present 2D radiative hydrodynamical simulationswhich show that multi-dimensional mixing of the super-nova ejecta with the nebular material may aid the es-cape of high-energy radiation. This could result in thegamma-ray escaping earlier than predicted by our calcu-lations, which assume a homogeneous ejecta shell. Po-larization studies of hydrogen-poor SLSNe at early timesfollowing the explosion show little evidence for signifi-cant asphericity of the outer ejecta layers (Leloudas et al.2015; Saito et al. 2020); however, higher polarizationshave been measured during the later nebular phase re-vealing more aspherical inner ejecta (Inserra et al. 2016;Saito et al. 2020). Nevertheless, the implied deviationsfrom spherical symmetry are not likely to alter our re-sults at the qualitative level.By contrast, several observations indicate that inAT2018cow the ejecta was highly aspherical, warrantingmore attention to this issue when expanding the tech-niques in this paper to the FBOT case. The opticaland X-ray light curves and spectra of AT2018cow canbe understood by at least two distinct ejecta compo-nents: a fast polar outflow of velocity ∼ > . ∼ < − (e.g. Margutti et al. 2019). In particu-lar, the lower densities and high ionization state of thefast polar outflow may have enabled soft X-rays to es-cape from the polar channel even at early times aroundoptical peak, while absorption of the same central X-raysource by the dense slower equatorial ejecta was respon-sible for the optical emission (see also Piro & Lu 2020). CONCLUSIONSWe have presented the first radiative transfer simu-lations of supernova light curves powered by the rota-tional wind of a pulsar/magnetar that self-consistentlycalculate the escape of high-energy radiation from thenebula and its thermalization by the ejecta, taking intoaccount a wide range of physical processes responsiblefor exchanging energy between the matter and radiationfields.Our results can be summarized as follows: • While some energy from the central engine maybe transferred directly to the ejecta at early timesafter the explosion (e.g., via mechanical work of ashock driven into the ejecta by the nebula; Kasenet al. 2016), the bulk of the optical radiation fromSLSNe around after peak light must be powered4
Vurm & Metzger by reprocessing of nonthermal radiation from thenebula by the supernova ejecta.Several processes in the ejecta are capable of ab-sorbing high-energy photons (Fig. 3), but an initialabsorption of energy is no guarantee of its ulti-mate thermalization. Thermalization requires re-processing the engine’s luminosity into low-energyphotons (which readily share their energy withthe thermal pool by photoelectric absorption orCompton down-scattering) or low-energy particles(which readily share their energy with the ejectaby Coulomb scattering). This behavior cannot becaptured by a single constant thermalization opac-ity for the nonthermal photons. • Gamma-rays are produced in the nebula by a com-bination of IC scattering and synchrotron emis-sion, which dominate the nebula luminosity atearly and late-times, respectively (the cross-overtime t ∼ t B ; eq. 17). Synchrotron emission tendsto produce lower energy photons, which are morereadily absorbed and thermalized by the ejecta,while IC emission produces higher energy gamma-rays which more readily escape without thermal-izing. The relative importance of these processesis sensitive to the nebula magnetization ε B , whichturns out to be the most important free variablein the problem (once the engine and ejecta prop-erties have been chosen to match the optical lightcurve near maximum). • Except in the case of very low magnetization ε B ∼ < − , the nebula remains radiative for yearsto decades or longer following the explosion, i.e.the bolometric output tracks the spin-down powerof the pulsar. However, the spectral energy dis-tribution of the nebular radiation depends sen-sitively on the mean per-particle energy ∆ γ (cid:39) L e / ( ˙ N ± m e c ) pairs gain upon entering the neb-ula from the pulsar wind. The latter depends onratio of the wind luminosity to the pair-loadingrate, ˙ N ± .Although the contribution to ˙ N ± that arises di-rectly from the inner magnetosphere is uncertain,this contribution is subdominant. Over the firstseveral years after the explosion, ˙ N ± is dominatedby pairs generated by γγ interactions in the up-stream wind region. We identify a new feedbackmechanism between the nebula radiation and γγ pair creation in the wind that regulates the meanenergy of the pairs entering the nebula. At veryearly and late times, nonthermal photons gener-ally serve as the targets for pair production, while the targets are thermal at intermediate times whenthe ejecta is becoming transparent to the highestenergy gamma-rays (Table 2, Appendix D). Thisself-regulation process has the benefit of renderingour results insensitive to the uncertain intrinsicpair-multiplicity of the wind.On the other hand, the reduced per-particle energyin the pulsar wind due to γγ mass-loading (oc-curring over timescales during which most of thespin-down power is released) may also negativelyimpact the efficacy of baryon cosmic ray acceler-ation and associated high-energy neutrino emis-sion in very young magnetar winds (e.g., Arons2003; Murase et al. 2015; Kotera et al. 2015; Fang& Metzger 2017; Fang et al. 2019). However, aquantitative exploration of the implication of ourresults for ion acceleration will require a more de-tailed model for the dissipation in the wind. • For the first several months after the explosion, in-cluding and following optical maximum, the ejectais opaque across all photon energies. During thisphase the engine’s luminosity is efficiently repro-cessed into thermal optical radiation, independentof any assumptions about the nebula. This pro-vides a rigorous justification for the assumptionof many previous works which calculate detailedoptical/UV light curves and spectra of magnetar-powered SLSNe around and following peak lightby depositing 100% of the magnetar’s spin-downluminosity as thermal energy at the inner edge ofthe ejecta (e.g. Kasen & Bildsten 2010; Dessartet al. 2012; Inserra et al. 2013; Metzger et al. 2015;Nicholl et al. 2017b). Furthermore, these modelscan indeed be used to obtain accurate constraintson the engine properties, such as the dipole mag-netic field strength and birth spin period, providedthat the engine lifetime t e (i.e. magnetar dipolespin-down time) is shorter than a few months. • Gamma-rays leak out of the ejecta starting monthsafter the explosion, beginning in the ∼ . − Fermi LAT and then mov-ing up to the ∼ >
100 GeV band accessible to at-mospheric Cherenkov telescopes. The shift in thegamma-ray spectra to higher energies is driven bya combination of the rising of mean particle energy∆ γ entering the nebula as the γγ pair-loading ofthe pulsar wind subsides and the opacity windowthrough the nebula and ejecta expands to higherphoton energies (both processes being driven bythe decreasing γγ optical depth, in the pulsar windand ejecta, respectively). amma-ray thermalization and leakage from millisecond magnetar nebulae ε B producing higher peak gamma-ray luminositiesand, conversely, lower late-time optical luminosi-ties. This is due to the tendency of more magne-tized nebulae to radiate a greater fraction of theengine’s energy in low-frequency synchrotron radi-ation, which is more readily absorbed by the ejectathan the higher frequency IC emission. • Reproducing the steep post-maximum decaysin the optical light curves of SN2015bn andSN2017egm requires weakly magnetized nebulae, ε B ∼ < − − − . These magnetizations aremuch lower than inferred from other young pul-sar wind nebula like the Crab Nebula. How-ever, the physical conditions in the extremelyyoung, rapidly spinning and highly magnetizedpulsars/magnetars considered here, may substan-tially differ from those of older pulsar nebulae fre-quently observed in the Milky Way. The partic-ular impact of γγ pair loading on potentially re-ducing the magnetization of the ultra-relativisticpulsar/magnetar wind feeding the nebula, meritsfurther study. • Alternatively, if such a low nebula magnetizationis deemed to be unphysical, our results suggestthat the magnetar model as currently envisionedmay be an incorrect or at least incomplete ex-planation for SLSNe. If the true luminosity ofthe central engine were to decrease faster in time(e.g., L sd ∝ t − . in SN2015bn) than the canonical ∝ t − magnetic dipole spin-down rate (e.g., due toa growing magnetic dipole moment or the effects offall-back accretion on the magnetar wind; Metzgeret al. 2018a), then even a sustained high level ofthermalization (as achieved in our high ε B models)would be consistent with the late-time optical lightcurve data. In such a model there would also beno heretofore unobserved large escaping gamma-ray flux.Absent such alternatives, other models such asCSM shock interaction or mildly-relativistic accre-tion disk outflows, could supplant ultra-relativisticmagnetar winds or black hole jets as the favoredengines of SLSNe. • Our models predict a late-time flattening in theoptical light curve, once synchrotron overtakes ICemission in the nebula on a timescale ∼ t B (eq. 17)and the effective thermalization opacity begins torise (Fig. 10). Hints of such flattening behaviorare evident in SN 2015bn (Fig. 6). • A definitive confirmation of the central enginemodel for SLSNe would come from the detectionof the leaking nebular gamma-rays. However, thisis an observational challenge due to the limitedsensitivity of gamma-ray telescopes and the de-layed onset of the gamma-ray flux, which reducetheir luminosity relative to the optical peak. Nev-ertheless, our models predict the expected gamma-ray emission given as the output of models whichcan self-consistently reproduce the observed opti-cal light curve. A detection may be possible for aparticularly nearby future SLSNe with
Fermi
LATat ∼ GeV energies, or present and future atmo-spheric Cherenkov telescopes such as the CTA at ∼ TeV energies. • Although our models are focused on SLSNe, thegeneral scenario we have outlined could havebroader applicability to other engine-poweredtransients such as FBOTs, tidal disruption events,and post-merger counterparts of neutron starmergers. All else being equal, the prospects maybe better for detecting the escaping high-energyemission from these events due to the shortertimescales over which the ejecta becomes transpar-ent to gamma-rays (Section 4.5). Indeed, the X-ray/gamma-ray emission from a magnetar or blackhole powered nebula may have already been de-tected in AT2018cow (e.g., Margutti et al. 2019).We thank Ke Fang, Keiichi Maeda, Matt Nicholl, andAkihiro Suzuki for providing helpful information andcomments on an early draft of the paper. IV acknowl-edges support by the ETAg grant PRG1006 and by EUthrough the ERDF CoE grant TK133. BDM acknowl-edges support from the NASA Astrophysics TheoryProgram (grant number NNX17AK43G), Fermi GuestInvestigator Program (grant number GG016287) andthe NSF through the AAG Program (grant numberGG016244).APPENDIX6
Vurm & Metzger A. RADIATIVE DIFFUSION WITH INCOHERENT SCATTERINGIn this Appendix we provide an approximate analytical prescription for computing the fraction of the energy injectedinto the ejecta by the central engine that emerges without significant attenuation. Neglecting emission from secondarypairs generated by the high-energy photons, the main remaining nontrivial issue in the problem is the characterizationof incoherent (Compton) scattering which degrades the photon energy without destroying it.To keep the problem simple while retaining the salient features, we treat Compton recoil as an energy sink rather thana frequency redistribution mechanism. All spectral information is lost in this approach, which however is acceptablefor the present purpose. This simplification allows us to consider radiative transfer independently at each frequency(without coupling between frequencies). The equation of radiative diffusion is given by dIds = − αI + j + (1 − κ C ) σJ, (A1)where I is the specific intensity, α is the absorption coefficient, j is the emissivity, σ is the scattering coefficient, J = (4 π ) − (cid:82) π Id Ω is the angle-averaged intensity, κ C is the fraction of the photon energy lost to recoil in a scatteringevent. Since we are considering high-energy radiation in an essentially cold medium, the emissivity j can be neglectedcompared with the last term in Equation (A1).Rewriting the transfer equation in spherical coordinates and taking the first two angular moments (where µ ≡ cos θ )yields 1 r ∂∂τ (cid:0) r H (cid:1) = − εJ (A2) ∂K∂τ − τ ( J − K ) = − H (A3)where H = (4 π ) − (cid:82) π Iµd
Ω and K = (4 π ) − (cid:82) π Iµ d Ω are the first and second moments of the intensity, respectively, τ = ( α + σ ) r , and ε = α + κ C σα + σ (A4)is defined such that 1 − ε is the effective single-scattering albedo. In this simplified treatment (cid:15) is the average fractionof the photon energy lost (either by absorption or recoil) in a single interaction.Following standard treatments (e.g. Rybicki & Lightman 1979) we use the Eddington approximation K ≈ J/ ∂ J∂τ (cid:63) + 2 τ (cid:63) ∂J∂τ (cid:63) − J = 0 , (A5)where τ (cid:63) = √ ε τ . Equation (A5) can be solved by standard methods, yielding J = C e τ (cid:63) τ (cid:63) + C e − τ (cid:63) τ (cid:63) . (A6)The integration constants C and C can be found by prescribing the flux at the inner boundary r in and requiringthat no radiation impinges on the outer boundary from above. The escaping flux can then be obtained from Equation(A3) (again using using J = K/ H = − (1 / ∂J/∂τ . The ratio of escaping energy to that input at the center of thespherical cloud is given by r H ( r ) r H ( r in ) | r in → = 2 τ (cid:63), [( √ (cid:15) − τ (cid:63), + √ (cid:15) ] e − τ (cid:63), + [( √ (cid:15) + 1) τ (cid:63), − √ (cid:15) ] e τ (cid:63), , (A7)where τ (cid:63), is the effective optical thickness of the cloud. Using the appropriate absorption and scattering opacities(see Section 2.3) in Equation (A7) yields the approximate fraction of energy that escapes without significant at-tenuation/reprocessing, for any given input photon energy. Figure 3 compares the effective opacity versus photonenergy derived from this result (dashed lines) compared to the result of a full Monte Carlo calculation (solid lines),demonstrating reasonable agreement at high values of τ T (cid:29) amma-ray thermalization and leakage from millisecond magnetar nebulae Figure 13.
Ratio of the dynamical timescale to cooling timescale in the supernova ejecta as a function of electron/positronLorentz factor γ , shown at three epochs after explosion ( t = 0.3 yr, 1 yr, 3 yr as marked). The dominant cooling process foreach range of γ is denoted by the line color as marked. We have assumed magnetar and ejecta parameters corresponding toSN2017egm (Table 3). The key point of this figure is to illustrate that relativistic electron/positrons are fast-cooling in theejecta across all relevant particle energies for several years after the explosion.B. PAIR COOLING TIMESCALESThe photon absorbing processes described in Section 2.3 result in the generation of energetic electrons or positrons,which cool and emit secondary radiation over a range of frequencies. The efficiency of this energy reprocessing dependson the ability of pairs to cool efficiently, i.e. faster than the expansion time; the relative importance of the differentcooling mechanisms depends on the radiation and matter field densities and has an impact on both the high-energyspectral evolution as well as the fraction of energy that ultimately gets thermalized. This Appendix covers the coolingrates of relativistic leptons of Lorentz factor γ due to different processes and compares the corresponding timescaleswith each other as well as the dynamical time t dyn = R ej /v ej over which adiabatic losses occur due to the expandingejecta/nebula.For relativistic bremsstrahlung emission, the particle cooling rate is ˙ γ br ≈ cα fs τ T γ . ( Z /Z + 1) / (6 R ej ) (Vurm &Metzger 2018). The ratio of the expansion to cooling timescales is t dyn t br = ˙ γ br γ R ej v ej = 5 cα fs τ T γ . v ej (cid:32) Z Z + 1 (cid:33) ≈ . γ . M v − t − (B8)where in the final step γ ≡ γ/ and we have assumed oxygen-rich ejecta, Z = ( Z ) / = 8.8 Vurm & Metzger
For IC cooling in the Thomson regime ( γ ∼ < ), for which the particle cooling rate is ˙ γ IC = 4 c(cid:96) th γ / (3 R ej ), weobtain t dyn t IC = ˙ γ IC γ R ej v ej = 4 c (cid:96) th γ v ej = 0 . γ M v − B − t − (B9)where (cid:96) th is the thermal compactness (eq. 10) and in the final step we have assumed κ γ /κ es = 0 . τ opt ≤ Finally, cooling due to Coulomb scattering ˙ γ Coul = 3 c ln Λ τ T / (2 R ej ), results in t dyn t Coul = ˙ γ Coul γ R ej v ej = 3 ln Λ c τ T v ej γ = 1 . γ − M v − t − , (B10)where in the final step we have taken ln Λ ≈
25 for the Coulomb logarithm. As long as the sum of (B8) - (B10) exceedsunity for a given set of parameters, the high-energy electron/positron is in the fast-cooling regime and its adiabaticlosses can be neglected in the lowest order.The relative rates of bremsstrahlung and Coulomb cooling can be expressed as˙ γ br ˙ γ Coul = (cid:18) γγ (cid:63) (cid:19) . , (B11)where γ (cid:63) = (cid:20) α fs ( Z /Z + 1) (cid:21) . ≈ (cid:32) Z + Z Z (cid:33) − . ≈
230 (B12)and in the last step we have again assumed oxygen-dominated composition.Similarly, for bremsstrahlung and IC one obtains˙ γ IC ˙ γ br = (cid:18) γγ (cid:63)(cid:63) (cid:19) . , (B13)where γ (cid:63)(cid:63) = (cid:34) τ T α fs ( Z /Z + 1)8 (cid:96) th (cid:35) . ≈ . × (cid:96) − . , − τ . (cid:32) Z + Z Z (cid:33) . ≈ . × (cid:96) − . , − τ . (B14)and (cid:96) th , − = (cid:96) th / − .Finally, the ratio of IC and Coulomb cooling rates is˙ γ IC ˙ γ Coul = (cid:18) γγ † (cid:19) , (B15)where γ † = (cid:18) τ T (cid:96) th (cid:19) / ≈ . × τ / l − / , − . (B16)Figure 13 shows the ratio of the dynamical timescale to the total cooling timescale in the supernova ejecta as a functionof γ at three epochs after explosion ( t = 0.3 yr, 1 yr, 3 yr) as calculated from the above estimates for engine/ejectaparameters corresponding to SN2017egm. At all epochs Coulomb losses dominate the energy loss of low- γ leptons, whilethe high- γ electrons cool through primarily IC emission at early times and via bremsstrahlung at later times. For allelectron energies γ adiabatic losses can be neglected for the first several years after the explosion (i.e. t dyn /t cool > The target radiation for IC cooling should in fact include allphotons in the Thomson regime, i.e. with x ∼ < /γ includingUV/X-rays. However, as those photons are efficiently photoelec-trically absorbed by the ejecta, the error made by using just thethermal radiation density is modest. amma-ray thermalization and leakage from millisecond magnetar nebulae C. RADIATIVE COOLING IN THE PULSAR WINDThis Appendix addresses whether the relativistically hot pairs deposited by γγ interactions in the pulsar/magnetarwind (Section 2.5) are able to cool radiatively before reaching the termination shock, or whether they will enter thenebula with their thermal energy intact and hence emit the engine luminosity isotropically as assumed in our MonteCarlo simulations (Section 3.1).Consider an electron or positron at the termination shock with an energy γ that emits IC photons on thermal targetradiation of temperature θ . Well after the optical peak the self-regulation mechanism of pair loading generally ensuresthat 4 θγ > x ≈ γ . A fraction of the upscatteredphotons propagates back into the unshocked wind and produces pairs via γγ interactions on the thermal radiation.The average energy of the generated electron-positron pair is γ ± , ≈ x/ ≈ γ /
2. If the pulsar wind is relativistic(bulk radial Lorentz factor Γ (cid:29)
1) and the generated electron/positron has time to isotropize in the wind frame, itsenergy in the lab frame after being picked up by the wind is γ ± ≈ Γ γ (cid:48)± ≈ Γ γ ± , ≈ Γ γ ≈ Γ x , (C17)where a prime denotes the wind rest frame.First consider IC cooling. The IC cooling time within the ambient thermal radiation field is (in the lab frame) t IC = 3 R t c(cid:96) th , w γ ± F KN , (C18)where (cid:96) th , w is the thermal compactness in the pulsar wind (see below) and F KN ≈ θγ ± ) (cid:20) ln(4 θγ ± ) − (cid:21) . (C19)is the Klein-Nishina correction factor for 4 θγ > t IC = R t c(cid:96) th , w x Γ (4 θγ ) ln(4 θγ ± ) − / . (C20)The cooling time (C20) should be compared with the time it takes for the generated pair to be advected back to theshock. On average, a high energy photon will propagate a distance ∆ = min( R t , R t /τ γγ, th ) back into the wind zonebefore producing a pair, so the average advection time is ∆ /c . The pair production opacity is τ γγ, th ≈ (cid:96) th , w x ln(1 + xθ )5 . xθ ) exp (cid:18) − . θx (cid:19) . (C21)Substituting (cid:96) th , w x from Equation (C21) back into Equation (C18), using Equation (C17) and rearranging terms, weobtain t IC t adv ≈ Γ ln(1 + θγ ) e − . /θγ ln(4 θγ ) + ln(Γ / − / , (C22)where t adv = ∆ c = R t cτ γγ, th (C23)and we have assumed that τ γγ, th > θγ >
1, the last fraction in Equation (C22) is comparable to unity, so we conclude that in a relativisticwind the generated pairs do not have time to cool by inverse Compton emission before being advected back to theshock.Now consider synchrotron cooling. The lab-frame synchrotron cooling time is t syn = Γ t (cid:48) syn = Γ 3 m e c σ T U (cid:48) B γ (cid:48)± = 3 m e c Γ σ T U B γ ± , = 3Γ (cid:96) B, w γ ± , R t c , (C24)0 Vurm & Metzger where again γ ± , is the energy of the created pair before being picked up by the wind. The ratio of cooling andadvection times is t syn t adv = 3Γ (cid:96) B, w γ ± , (cid:96) th , w x ln(1 + xθ )5 . xθ ) exp (cid:18) − . θx (cid:19) ≈ (cid:96) th , w (cid:96) B, w Γ ln(1 + θγ ) e − . /θγ (4 θγ ) , (C25)where again τ γγ, th > (cid:96) B, w = σ T U B R t m e c = (cid:96) inj σσ + 1 R ej R t , (C26)where U B = B / π , and the magnetic field and injection compactness are (eq. 20) B w = (cid:20) L e σ ( σ + 1) cR (cid:21) / (C27)and (cid:96) inj = σ T L e πm e c R ej , (C28)respectively. The thermal compactness in the wind zone is (cid:96) th , w = σ T R t u th , w m e c ≈ σ T R t m e c L e f th (1 + τ opt )4 πcR = (cid:96) inj f th (1 + τ opt ) R ej R t . (C29)Using the above definitions, the ratio of the synchrotron cooling and advection times becomes (again for τ γγ, th > t syn t adv ≈ σ + 1 σ f th (1 + τ opt )Γ ln(1 + θγ ) e − . /θγ (4 θγ ) . (C30)Given that pair creation in the wind regulates θγ ∼ few (Section 2.5), we see that even for a strongly magnetizedwind σ (cid:29)
1, synchrotron cooling is negligible ( t syn (cid:29) t adv ) provided that the wind is relativistic Γ ∼ > f − / ∼ − τ γγ, th < t adv = R t /c , the critical ratio instead becomes t syn t adv = 3Γ (cid:96) B, w γ ± , = 3Γ (cid:96) inj γ σ + 1 σ R t R ej , (C31)which, by definition, must be smaller than the value given by Equation (C30). In this case, synchrotron cooling is onlynegligible if the wind remains highly relativistic and/or the wind magnetization has dropped to low values by the timeit reaches the pair loading zone. However, once τ γγ, th < γγ pair loading on the radial evolution of the Lorentz factor andmagnetization of the pulsar wind (see §
4) and to assess the conditions required for particle cooling as described inthis section within such a self-consistent framework. However, the results derived here make it plausible that radiativelosses in the wind may be neglected to first order, supporting the assumptions of our simulations (Section 3.1). D. REGULATION OF THE PARTICLE NUMBER AND MEAN PARTICLE ENERGY IN THE WINDDISSIPATION REGIONUpon entering the nebula, a heated lepton cools rapidly by both IC scattering on thermal photons and synchrotronemission (Section 2.2), with a fraction ξ ≈ / − e − τ γγ . Putting this together withEquation (34) yields the following condition: M e = 2 ξ (cid:90) dN ph d ln x (cid:12)(cid:12)(cid:12)(cid:12) . (1 − e − τ γγ ) d ln x = 1 , (D32) amma-ray thermalization and leakage from millisecond magnetar nebulae dN ph /d ln x is thespectrum emitted by a single electron over its cooling history.To make further progress on must specify the spectrum emitted by the electron as it cools down from its initial γ ≈ ∆ γ , as well as the source of γγ opacity. It is useful to discuss separately the cases when the opacity is providedby nonthermal or thermal target photons, as well as to separate the cases based to the physical mechanisms givingrise to the high-energy and target photons. Results for these different cases are summarized in Table 2.D.1. Nonthermal pair loading
D.1.1.
IC photons + IC targets
Close to the optical peak the injection compactness (cid:96) inj ∼ . − ∼ < γ (Vurm & Metzger 2018) dN ph d ln x dt (cid:12)(cid:12)(cid:12)(cid:12) IC , . = ˙ γ IC δ (cid:0) x − θγ (cid:1) . (D33)The spectrum emitted by the electron as it cools down from its initial γ ≈ ∆ γ is obtained by integrating Equation(D33) over time dt = dγ/ ˙ γ , yielding dN ph d ln x (cid:12)(cid:12)(cid:12)(cid:12) IC , . = 12(4 θx ) / , x ≤ θγ , (D34)where we have assumed ˙ γ ≈ ˙ γ IC . The deposited (but not necessarily escaping) IC luminosity per logarithmic photoninterval is dLd ln x = ˙ N ± x dN ph d ln x (cid:12)(cid:12)(cid:12)(cid:12) IC , . = L e xγ dN ph d ln x (cid:12)(cid:12)(cid:12)(cid:12) IC , . , (D35)where ˙ N ± denotes the number of pairs, including secondaries, that share the dissipated power. The frequency-dependent compactness (eq. 29) that enters opacity (eq. 28) can now be written as (cid:96) ν = σ T m e c R n du γ d ln x = σ T m e c R n dLd ln x t res V n = 3 σ T m e c L e πR n xγ dN ph d ln x (cid:12)(cid:12)(cid:12)(cid:12) IC , . t res t LC = 32 (cid:96) inj R ej R n (cid:18) xx IC (cid:19) / t res t LC , (D36)where t LC = R n /c is the light crossing time of the wind/nebula, x IC = 4 θγ is the characteristic IC photon energy, andthe injection compactness (cid:96) inj is defined by equation (8). The target photon residence time t res in the wind/nebula istypically of the same order as t LC , since the nebula is baryon starved and both photoionization and Compton recoillosses are negligible.Using equations (28), (D34) and (D36) in equation (D32) and assuming τ γγ, nth < M e ≈ ξ η ( α ) (cid:96) inj R ej R n ln x IC θγ = 1 . (D37)A few remarks should be made about equation (D37). First, equation (D37) can only be satisfied if x IC = 4 θγ > x IC ∼ >
1, the2
Vurm & Metzger term (4 θγ ) − = (4 θx IC ) − / = 500 θ − / − x − / (cid:29)
1, so that M e > γ > (4 θ ) − / as long as (cid:96) inj remains sufficiently high. Therefore, to satisfy M e = 1, the pair loading maintains γ sufficiently low so that only asmall fraction of IC photons exceed the pair production threshold. In other words, x IC ≈ γ ≈ (4 θ ) − / , as longas (cid:96) inj ∼ > (4 θ ) / η ( α ) R n R ej ≈ . θ / − R n R ej , (D38)where we have used η ( α ) ≈ .
1. Using equation (8) for (cid:96) inj , this occurs on a timescale t ± , IC ≈ . B − / v − / θ − / − yr . (D39)Finally, once (D38) is no longer satisfied ( t > t γγ, nth ), the nonthermal pair regulation rapidly breaks down: since M e ∝ γ − , the feedback process is unstable away from the threshold (at x IC (cid:29) N ± produces a small increase in γ , which leads to a decrease in M e and hence a further decrease in ˙ N ± . The feedbackthen switches to the thermal mechanism described in the next subsection.D.1.2. Synchrotron photons + synchrotron targets
Another mode of nonthermal pair loading can occur if at some stage synchrotron photons become sufficiently energeticfor pair creation, i.e. if the dimensionless synchrotron frequency x syn = hν syn , / ( m e c ) = γ x B ≥
1, where x B is thedimensionless Larmor frequency. The argument is analogous to the one presented above, except one has to replace˙ γ IC , 4 θ and x IC with ˙ γ syn , x B and x syn , respectively, in equations (D33)-(D37). This yields a condition analogous to(D37), under the assumption that pair cooling is dominated by sychrotron radiation. One can account for additionalIC losses on thermal radiation by writing ˙ γ = ˙ γ syn + ˙ γ IC = ˙ γ syn (1 + (cid:96) th F KN /(cid:96) B ), where F KN ≈ , θγ (cid:28) θγ ) (cid:20) ln(4 θγ ) − (cid:21) , θγ (cid:29) . (D40)is the Klein-Nishina correction to the Thomson IC cooling rate. Using this in dt = dγ/ ˙ γ when integrating equation(D33) (with the above substitutions) over the electron cooling history, the synchrotron spectrum emitted by a singleelectron becomes dN ph d ln x (cid:12)(cid:12)(cid:12)(cid:12) syn , . = 12( x B x ) / (1 + (cid:96) th F KN /(cid:96) B ) , x ≤ x B γ . (D41)Instead of condition (D37) one now obtains M e ≈ ξ η ( α ) (cid:96) inj (cid:96) th F KN /(cid:96) B R ej R n ln x syn x B γ = 1 . (D42)The threshold condition hν syn , > m e c requires (eq. 26) γ > γ † = 6 . × B / ε − / , − v / t yr . (D43)Below we will show that thermal pair loading regulates the injected lepton energy to γ ≈ /θ ≈ , hence thesynchrotron mechanism is unlikely to operate as long as thermal pair loading is effective. Once the thermal regulationfails, the synchrotron loading can operate as long as the “unloaded” γ ± = γ in (eq. 18) satisfies the above thresholdcondition.Neglecting IC cooling, the pair multiplicity (D42) can be written approximately as M e ∼ . (cid:96) inj γ ln x syn x syn = 0 . (cid:96) inj γ † γ † ln( γ /γ † ) γ , γ ≥ γ † , (D44)where we have used γ † x B = 1. Intuitively, the multiplicity M e is a product of γ /x syn ∼ γ photons produced perlepton and a fraction τ γγ ∼ . (cid:96) inj of them producing secondary pairs. The second equality in (D44) shows that as amma-ray thermalization and leakage from millisecond magnetar nebulae . (cid:96) inj γ † > M e ( γ ) = 1 has two roots, γ ∼ γ † and another at γ > γ † .The latter solution is in the unstable region, where M e ∝ γ − , i.e. any initial perturbation of γ is amplified untilanother mechanism takes over the pair regulation. Typically, γ is a growing function of time, so the lower-energy(stable) solution γ ∼ γ † is attained first; the system remains in that state as long as 0 . (cid:96) inj γ † > t ± , syn obeying t ± , syn < min( t (2) ± , syn , t (1) ± , syn ) . (D45)The first timescale is the requirement that 0 . (cid:96) inj γ † >
1, which is satisfied before a time t (1) ± , syn ≈ . B − / v − / ε − / B, − yr . (D46)The second timescale arises from the condition that the regulated γ ∼ γ † be lower than the “unregulated” value γ in (eq. 18), which yields t (2) ± , syn = 7 . B − / v − / ε / B, − (cid:16) µ ± (cid:17) − / yr . (D47)D.1.3. IC photons + synchroton targets
The mechanisms outlined in the previous two sections operate in the early (close to the optical peak) and late(years) stages of the event, respectively. In addition to the thermal mechanism discussed below, a further mode of pairregulation is possible in the intermediate stage ( t ∼ x IC x syn = 4 θx B γ >
1. The multiplicity M e isobtained from equation (D32) by using equation (D34) for the IC spectrum emitted by a single electron (with anadditional factor (1 + (cid:96) B /(cid:96) th F KN ) − that accounts for synchrotron losses), and using the synchrotron spectrum (D41)for computing the pair production opacity. The result is M e ≈ ξ η ( α ) (cid:96) inj (cid:96) B (cid:96) th F KN ( (cid:96) B + (cid:96) th F KN ) R ej R n γ ln( x IC x syn )( x syn x B ) / = 1 . (D48)Following an argument analogous to the previous two sections, one concludes that pair loading maintains γ ∼ (4 θx B ) − / (i.e. x IC x syn ∼ γ ≈ . × θ − / − B / ε − / , − v / t / , (D49)as long as 0 . γ (cid:96) inj (cid:96) B (cid:96) th / ( (cid:96) B + (cid:96) th ) > F KN ≈ t < t ± , IC − syn , where t ± , IC − syn ≈ θ − / − B − / ε − / , − v − / (cid:20) (cid:96) B (cid:96) th ( (cid:96) B + (cid:96) th ) (cid:21) / yr . (D50)D.2. Thermal pair loading
In this regime pair loading is controlled by IC gamma-rays interacting with thermal target photons. To determinethe pair multiplicity, we now need to use the IC spectrum of an electron cooling on the thermal radiation field in bothThomson and Klein-Nishina regimes. In the delta-function approximation, one can write dN ph d ln x dt (cid:12)(cid:12)(cid:12)(cid:12) . = (cid:40) ˙ γ IC δ (cid:0) x − θγ (cid:1) , θγ ≤ γ IC δ ( x − γ ) , θγ > . (D51)The total cooling rate, also including synchrotron emission, can be written as ˙ γ = ˙ γ IC + ˙ γ syn = ˙ γ IC (cid:18) (cid:96) B (cid:96) th F KN (cid:19) = 4 c R n (cid:96) th γ F KN (cid:18) (cid:96) B (cid:96) th F KN (cid:19) . (D52) We have neglected Coulomb and bremsstrahlung cooling (Ap-pendix B), as justified by the low Thomson optical depth andlow expected baryon content of the wind and nebula. Vurm & Metzger
The spectrum emitted by a cooling electron is obtained by integrating Equation (D51) over time dt = dγ/ ˙ γ : dN ph d ln x (cid:12)(cid:12)(cid:12)(cid:12) . = θx ) / (1 + (cid:96) B /(cid:96) th ) , x ≤ min(1 / θ, θγ )11 + (cid:96) B / ( (cid:96) th F KN ) , / θ < x < γ , else. (D53)The pair production opacity τ γγ, th that enters Equation (D32) drops rapidly for 4 θx < θγ ∼ >
1, i.e. the cooling electron initially scatters thermal photons in theKlein-Nishina regime. On the other hand, the number of emitted IC photons increases with decreasing photon energy.Therefore, the main contribution to the integral in eq. (D32) comes from values of 4 θx ≈ few.Using the spectrum (D53) we can rewrite condition (D32) as: M e ≈ ξ (cid:90) / θ θx ) / (1 + (cid:96) B /(cid:96) th ) (1 − e − τ γγ, th ) d ln x + 2 ξ (cid:90) γ / θ
11 + (cid:96) B / ( (cid:96) th F KN ) (1 − e − τ γγ, th ) d ln x = 1 . (D54)Figure 5 shows the numerically computed multiplicity M e , using now exact IC/synchrotron cooling rates and spectra,as a function of γ for different values of the gamma-ray optical depth τ γγ, th through the wind region (evaluated atthe peak value xθ ≈
2) and different ratios of the magnetic to thermal compactness, (cid:96) B /(cid:96) th (eq. 15). As expected, thepair multiplicity increases with the initial Lorentz factor of the radiating pair γ , the temperature of the backgroundradiation field θ , and the gamma-ray optical depth τ γγ, th . However, to understand the detailed behavior of M e withincreasing γ θ we must rely on the approximate expression (D54). Through condition (34) this determines the regulatedLorentz factor of the injected nebular pairs (eq. 21) and how it evolves as τ γγ, th decreases and (cid:96) B /(cid:96) th increases as timeadvances after the explosion.First, consider the limit of low (cid:96) B /(cid:96) th (cid:28) τ γγ, th >
1, which characterizes the nebula/wind region at earlytimes of months to years (depending on the value of ε B ; eq. 16). In this limit, the factor 1 − e − τ γγ, th ≈ xθ < τ γγ, th on xθ . As a result, the firstterm on the right hand side of Equation (D54) is typically less than unity unless the thermal compactness is very high.The second integral in Equation (D54) is trivial in the domain of interest as long as τ γγ, th (cid:29) (cid:96) B /(cid:96) th (cid:28)
1; thecondition M e = 1 yields ln(4 γ θ ) ∼
1, or 4 γ θ ∼ a few. Consistent with this, Fig. 5 shows that the value γ whichsatisfies M e = 1 (eq. 34) depends only weakly on τ γγ, th for τ γγ, th > (cid:96) B /(cid:96) th the integrands of both terms in Equation (D54) are smaller, such that a higher valueof γ is required to attain M e = 1. However, once (cid:96) B /(cid:96) th increases to a value ∼ > .
1, the pair multiplicity dropssignificantly and eventually the thermal regulation condition M e = 1 can no longer be satisfied for any γ . Physically,this is because synchrotron emission (which dominates IC cooling for (cid:96) B /(cid:96) th ∼ > .