A Generalized Semi-Analytic Model for Magnetar-Driven Supernovae

Several types of energetic supernovae, such as superluminous supernovae (SLSNe) and broad-line Ic supernovae (Ic-BL SNe), could be powered by the spin-down of a rapidly rotating magnetar. Currently, most models used to infer the parameters for potential magnetar-driven supernovae make several unsuitable assumptions that likely bias the estimated parameters. In this work, we present a new model for magnetar-driven supernovae that relaxes several of these assumptions and an inference workflow that enables accurate estimation of parameters from lightcurves of magnetar-driven supernovae. In particular, in this model, we include the dynamical evolution of the ejecta, coupling it to the energy injected by the magnetar itself while also allowing for non-dipole spin down. We show that the model can reproduce SLSN and Ic-BL SN light curves consistent with the parameter space from computationally expensive numerical models. We also show the results of parameter inference on four well-known example supernovae, demonstrating the model's effectiveness at capturing the considerable diversity in magnetar-driven supernova lightcurves. The model fits each light curve well and recovers parameters broadly consistent with previous works. This model will allow us to explore the full diversity of magnetar-driven supernovae under one theoretical framework, more accurately characterize these supernovae from only photometric data, and make more accurate predictions of future multiwavelength emission to test the magnetar-driven scenario better.

Several different models can be used to explain the high energies of one or both of SLSNe or SNe Ic-BL.Stars with masses M * ≳ 130 M ⊙ can explode as pair instability supernovae (PISNe) (Barkat et al. 1967;Heger & Woosley 2002;Gal-Yam et al. 2009), which can generate tens of solar masses of 56 Ni and present as an extremely long-lived, luminous supernovae.Slightly less massive stars, with 130 M ⊙ ≳ M * ≳ 100 M ⊙ , can eject shells of material through the pair-instability without being completely destabilized, and collisions between the ejected shells can be as luminous as an SLSN (Heger et al. 2003;Woosley et al. 2007;Chatzopoulos & Wheeler 2012;Yoshida et al. 2016;Woosley 2017) -these are known as pulsational pair-instability supernovae (PPISNe).The collision of supernova ejecta and circumstellar material surrounding the progenitor, ejected through either a steady wind, eruptive mass loss, or binary interaction (Smith 2014), can convert much of the supernova kinetic energy into radiation, leading to a highly luminous supernova (Chatzopoulos & Wheeler 2012; Chatzopoulos et al. 2013;Villar et al. 2017;Jiang et al. 2020).Finally, the compact remnant can inject some energy into the ejecta; this energy may come from fallback accretion onto a central black hole or neutron star (Dexter & Kasen 2013;Moriya et al. 2018), a collapsar or jet (MacFadyen & Woosley 1999), or the spin-down energy of a rapidly-rotating magnetar (Kasen & Bildsten 2010;Woosley 2010).
In the magnetar model, the spin-down energy is emitted as a highly magnetized particle wind, which expands relativistically until it collides with the inner edge of the supernova ejecta, create forward and reverse shocks.The expanding wind becomes shocked by the reverse shock, accelerating the particles up to ultrarelativistic energies, which then emit via synchrotron radiation and inverse Compton scattering (Gaensler & Slane 2006); this shocked wind is known as a pulsar wind nebula (PWN).The PWN applies a pressure to the supernova ejecta, causing it to accelerate, and the radiation from the nebula is thermalized in the ejecta, causing the ejecta temperature and supernova luminosity to both increase (Kasen & Bildsten 2010).The PWN-ejecta interaction can also cause Rayleigh-Taylor instabilities, which can shred the inner ejecta and cause non-spherical structure to emerge in the ejecta (Suzuki & Maeda 2017, 2021).
The magnetar model predicts many multiwavelength signals that could be used to identify and characterize the newborn neutron star.Once the ejecta becomes optically thin, the non-thermal emission from the PWN can be detected directly, either at high energy in hard X-rays or gamma rays (Kotera et al. 2013;Murase et al. 2015;Kashiyama et al. 2016), or at low energy in radio (Omand et al. 2018;Eftekhari et al. 2021;Murase et al. 2021); two energetic SNe have radio detections at late times that are consistent with PWN emission, PTF10hgi (Eftekhari et al. 2019;Law et al. 2019;Mondal et al. 2020;Hatsukade et al. 2021) and SN2012au (Stroh et al. 2021).Dust formed in the supernova can absorb PWN emission and re-emit that energy in infrared (Omand et al. 2019), causing bright continuum emission recently seen in four SLSNe (Chen et al. 2021;Sun et al. 2022).High-energy photons and Rayleigh-Taylor induced mixing can change the chemical and ionization structure of the ejecta, leading to unique signatures in the supernova nebular spectrum (Chevalier & Fransson 1992;Omand & Jerkstrand 2023).Aspherical ejecta caused by either hydro instabilities or an aspherical PWN can produce an optical polarization signal (Tanaka et al. 2017), which has been detected in some energetic SNe (Inserra et al. 2016;Saito et al. 2020;Pursiainen et al. 2022;Poidevin et al. 2023;Pursiainen et al. 2023), but not others (Leloudas et al. 2015b;Lee 2019Lee , 2020;;Poidevin et al. 2022Poidevin et al. , 2023;;Pursiainen et al. 2023).
Accurate parameter estimation from the light curve around optical peak is essential for a number of reasons.Firstly, new surveys such as the LSST will detect SLSNe out to high redshift, where they can not all be classified with spectroscopy.Therefore, these supernovae will have to be characterized from their light curve data alone.Also, predicting late-time multiwavelength signals requires accurate parameter estimation, as different parameters that produce similar optical light curves can produce vastly different multiwavelength signals (e.g.Omand et al. 2018Omand et al. , 2019)).In most widely used parameter inference codes (e.g Nicholl et al. 2017b), the model used for inference of magnetar-driven SNe makes several assumptions to reduce computational complexity, which are unjustified outside a small region of the parameter space.In particular, they assume a constant ejecta velocity, which is independent of the magnetar luminosity, although numerical simulations show that ejecta acceleration due to PWN pressure plays a vital role in the dynamics of the ejecta (Chen et al. 2016;Suzuki & Maeda 2017, 2019;Chen et al. 2020;Suzuki & Maeda 2021).They also assume the magnetar spin down through pure vacuum dipole emission, even though studies of Galac-tic pulsars (Lyne et al. 2015;Parthasarathy et al. 2020) and putative magnetars born in GRBs (Lasky et al. 2017;Sarin et al. 2020a,b) show that most neutron stars are inconsistent with a pure vacuum dipole.
In this work, we present a model where these assumptions are relaxed, which can fully explore the diversity of magnetar-driven supernovae and unite phenomenologically different supernovae, such as SNe Ic-BL and SLSNe, under one theoretical framework.In Section 2, we introduce our model for magnetar-driven supernovae.In Section 3, we show the diversity of supernovae and supernova observables resulting from differences in initial parameters.In Section 4, we perform Bayesian inference on a few varying supernovae to show how our model can consistently reproduce and explain them.Finally, in Section 5, we discuss the model's implications and conclude.Throughout the paper, we use the notation   = /10  in cgs units unless otherwise noted.

MODEL
The physics of our model is based on previous magnetar-driven kilonovae models (e.g.Yu et al. 2013;Metzger 2019;Sarin et al. 2022), but modified to describe supernovae.We present here a nonrelativistic model description, although the model implementation is fully relativistic.For the fully relativistic description of the kinematics (Equations 5,6,14,16,17,18,and 19), see Sarin et al. (2022).We note that relativistic corrections are largely unimportant for supernovae apart from transients with exceptionally low ejecta masses and powerful magnetar engines.

Model Physics
The central magnetar spins down by releasing its rotational energy where  is the moment of inertia of the magnetar and Ω is the rotational angular frequency of the magnetar.The time derivative of this relation gives the spin down luminosity, which, given Ω ∝ −Ω  for braking index , can be modelled generally as (Lasky et al. 2017) where  0 is the initial magnetar spin-down luminosity and  SD is the magnetar spin-down time.The vacuum dipole spin-down mechanism (Ostriker & Gunn 1969;Goldreich & Julian 1969) has a corresponding braking index of  = 3, which gives a late-time spin-down down luminosity of  ∝  −2 (Zhang & Mészáros 2001), while gravitational wave spin down via magnetic deformation (Cutler & Jones 2000) has a corresponding braking index of  = 5, and has a latetime spin-down luminosity of  ∝  −3/2 .There are also several other mechanisms for spin down, including particle winds (Harding et al. 1999;Xu & Qiao 2001;Wu et al. 2003;Contopoulos & Spitkovsky 2006), multipolar radiation (Pétri 2015), and r-mode oscillations (Ho & Lai 2000), among others (Blandford & Romani 1988;Lin & Zhang 2004;Ruderman 2005;Chen & Li 2006;Yue et al. 2007;Contopoulos 2007), and several mechanisms can spin-down the magnetar simultaneously.We note that in general,  is also expected to be variable during the early life of the magnetar (Lander & Jones 2018, 2020), but we keep it constant for simplicity.Integrating Equation 3 gives the total rotational energy as a function of braking index: This recovers  rot =  0  SD for vacuum dipole spin-down.The rotational energy from the magnetar is converted into a pulsar wind.This highly magnetized, ultrarelativistic wind collides with and pushes a shock into the supernova ejecta, increasing its kinetic and internal energy.The internal energy is also increased by the absorption of PWN photons by the ejecta.The total energy of the system is the combination of the kinetic and internal energy where  ej is the ejecta mass.The evolution of this system is governed by the energy sources, radioactive heating and magnetar spin-down luminosity, and energy losses from radiated luminosity and adiabatic cooling from expansion.The evolution of the internal energy of the ejecta is written as (Kasen et al. 2016) where  ra and  bol are the radioactive power and emitted bolometric luminosity, respectively,  is the fraction of spin-down luminosity injected into the ejecta, and P and  are the pressure and volume of the ejecta.
Here, we adopt the Wang et al. (2015) prescription for gamma-ray leakage used in other models (e.g.Nicholl et al. 2017b;Sarin et al. 2022) with where is the leakage parameter and   is the gamma-ray opacity of the ejecta.
The radioactive power from the decay of 56 Ni is given by where  Ni is the nickel fraction of the ejecta, 56 Ni = 6.45 × 10 43 erg s −1  −1 ⊙ and 56 Co = 1.45 × 10 43 erg s −1  −1 ⊙ are the decay luminosities of 56 Ni and 56 Co, and 56 Ni = 8.8 days and 56 Co = 111.3days are the decay timescales for 56 Ni and 56 Co (Nadyozhin 1994).The inclusion of  ra here allows our model to simplify to an Arnett (1982) model for low  SD .
The dynamical evolution of the ejecta is given by (Sarin et al. 2022) where Substituting these into Equation 10gives The initial ejecta velocity is set by where  SN is the supernova explosion energy.
The bolometric radiated luminosity is (Kasen & Bildsten 2010;Kotera et al. 2013) where is the optical depth of the ejecta,  is the ejecta opacity, is the effective diffusion time, and   >  dif is the time when  = 1.
Calculating the bolometric luminosity of the magnetar-driven transient (Equations 16 and 17) involves solving the evolution of the internal energy and dynamics (Equations 6 and 14 respectively) using the input power sources (Equations 3 and 9).The photospheric temperature is determined from the bolometric luminosity and ejecta radius until the temperature reaches the photospheric plateau temperature, as in Nicholl et al. (2017b).This can be expressed as The spectral energy distribution (SED) of the transient is then calculated using the cutoff blackbody used in Nicholl et al. (2017b), with  < cut =   (/ cut ) and  cut = 3000 (Chomiuk et al. 2011;Nicholl et al. 2017a), but this can also be switched to a simple blackbody.

Parameters, Priors, and Implementation
The model presented above is implemented into the open-source electromagnetic transient fitting software package, Redback (Sarin et al. 2023).The input parameters for the model, their default priors, and the values used in Section 3 are listed in Table 1.
The range of the default priors can lead to magnetars with rotational energies greater than 10 58 erg, which is unphysically high.When performing inference, the prior can be restricted to only sample parameters where  rot is below a certain threshold.We recommend using this constraint to prevent unphysically high energies, and use a value of 10 53 erg for the case studies presented in Section 4.
Due to the acceleration of the ejecta from the PWN,  ej is not constant and can not be used as a free parameter, since it is coupled to both  SN and  SD .Velocity information from spectroscopy can be used to weight the results, although the velocity measured from  3 × 10 3 , 10 4 ]  5000   Table 1.The parameters and default priors for the generalized magnetar model, as well as the values used for the parameter exploration in Section 3. Priors are either uniform (U) or log-uniform (L).
absorption widths is not the same as the photospheric velocity which is not the same as the ejecta velocity (Arnett 1982), so we caution against this unless the velocities are well calibrated (e.g.Dessart et al. 2016).
Our model differs from others (e.g.Nicholl et al. 2017b) in the choice of input parameters used to determine the magnetar luminosity.Previous models use the initial magnetar spin period  0 , dipole component of the pulsar magnetic field , and neutron star mass  NS ; while we use  0 and  SD (and , which is implicitly fixed to 3 in other models).The interpretation of braking index can be complicated due to the number of possible spin-down mechanisms, so users wanting to test specific mechanisms can swap different magnetar implementations due to the modular nature of the model and Redback.We note that as the magnetar luminosity for vacuum dipole spin-down can determined from only  0 and  SD (see Equation 3), using three parameters is unnecessary for parameter inference.Using these parameters also avoids assumptions such as the moment of inertia of the neutron star, which depends on the equation of state (EoS) (Lattimer & Schutz 2005) and can vary depending on the mass and spin period of the neutron star (Worley et al. 2008).To recover parameters such as the magnetar spin period and magnetic field, one can use the scalings used in other models such as for transients consistent with  = 3, which assumes a neutron star with moment of inertia of ∼ 1.3 × 10 45 g cm 2 for a 1.4  ⊙ neutron star, which is consistent with the APR EoS (Akmal et al. 1998) or MDI EoS with density independent nuclear symmetry energy (Das et al. 2003;Shetty et al. 2007), which both give neutron star radii ∼ 11.5 − 12 km (Worley et al. 2008).These scalings become more complicated for  ≠ 3 (e.g., Shapiro & Teukolsky 1983), with other dependencies such as the ellipticity of the neutron star or bulk viscosity (depending on the spin-down processes involved), leaving it difficult to definitively recover a spin period or magnetic field with such simplified models.

RESULTS
We now explore the diversity of magnetar driven-supernovae for different initial conditions using the model derived in Section 2. We assume the supernova explosion energy is 10 51 erg, typical for a neutrino-driven explosion, and the supernova is entirely powered by magnetar energy after the initial explosion, with no contribution from 56 Ni.We fix both the ejecta optical opacity and gamma-ray opacity to 0.1 cm 2 g −1 .This optical opacity value is typical of strippedenvelope supernovae (Inserra et al. 2013;Kleiser & Kasen 2014) and the gamma-ray opacity will not affect the peak light curve properties unless it is extremely low, so we use the same value for both opacities here for simplicity.We also fix the photospheric plateau temperature to 5000 K, which is typical for SLSNe (Nicholl et al. 2017b) and SNe Ic-BL (Taddia et al. 2019), although slightly lower than the temperature of 6000 K suggested by Nicholl et al. (2017b), which corresponds to the recombination temperature of O II; ultimately, this temperature discrepancy will have no impact on the peak light-curve properties in most cases, and will only affect the light curve at late times.
First, we show the diversity of supernovae that can be produced using this model in Section 3.1.To compare with previously derived results, we assume vacuum dipole spin-down and use Equations 21 and 22 with a 1.4  ⊙ neutron star to express our initial conditions in terms of  0 and ; the mapping between ( 0 ,  SD ) and ( 0 , ) is shown in Figure 1.Then we show the effect of changing braking index in Section 3.2.
Several inefficiencies can prevent all the magnetar spin-down luminosity from eventually being emitted as supernova luminosity.Some fraction ( from Equation 6) will escape without interacting with the ejecta at all, and possibly be detected as non-thermal x-rays or gamma rays.Some fraction of the energy will also accelerate the ejecta instead of thermalizing and being re-emitted, which will affect both the supernova luminosity and peak timescale; this fraction can be determined by the ratio of the spin-down time and the diffusion time (Suzuki & Maeda 2021;Sarin et al. 2022).
Figure 2 shows the ratio of the final supernova kinetic and radiated energies for various magnetic fields and ejecta masses for spin periods of 1 ms (close to the mass shedding limit for neutron stars (Watts et al. 2016)) and 3 ms (where the spin-down luminosity and explosion energy can become comparable), and for various magnetic fields and spin periods for an ejecta mass of 10  ⊙ , the same ejecta mass as Suzuki & Maeda (2021).The energy ratio  kin / rad does correlate strongly with  =  SD / dif up to large ejecta mass for low spin periods, although as spin period increases, the correlation gets weaker as the behaviour of the energy ratio changes; this is due to the total amount of energy injected by the magnetar prior to the diffusion time decreasing below the explosion energy, meaning that the ejecta dynamics are no longer primarily determined by the magnetar spindown energy.There is a small region of the parameter space at low spin period, magnetic field, and ejecta mass where the radiated energy can surpass the kinetic energy.However, for most of the parameter space this ratio is between 1 and 100.Typical supernovae have  kin / rad = 10 51 erg / 10 49 erg = 100, but without a contribution from 56 Ni this ratio can go much higher.Suzuki & Maeda (2021) find  kin / rad to be lower than we do in their 2D simulations by a factor of ∼ 2-10 depending on the model.However, they also find that multi-dimensional effects cause the supernova to be more luminous compared to 1D models at early times due to hot bubble breakout in the ejecta.
Figure 3 shows the peak timescale of the bolometric luminosity and the final ejecta velocity over the same parameter grid as Figure 2. The trends match our theoretical expectations, as  0 increases, the total energy of the system decreases, causing the ejecta velocity to decrease and the peak timescale to increase.Conversely, as  increases, the spin-down time decreases, causing the ejecta velocity to increase and the peak timescale to decrease.Finally, as  ej increases, the diffusion time increases, causing the ejecta velocity to decrease and the peak timescale to increase.The final ejecta velocity is largely insensitive to the magnetic field strength above a certain field threshold despite a change in peak timescale.This is a product of the coupling of dynamical evolution of the ejecta to the magnetar's rotational energy, as the final velocity is reached earlier in the supernova evolution for higher magnetic field (due to their shorter spin-down timescale).For a magnetar with rotation period of 1 ms, a magnetic field of > 10 15 is needed to reduce the peak timescale below 10 days, even at an ejecta mass of 1  ⊙ , meaning that the fastest SNe Ic-BL likely require both an ejecta mass below 1  ⊙ and a magnetar spinning at close to breakup speeds.Timescales of > 20 days require ejecta masses below 5  ⊙ , meaning that our model will likely estimate a lower ejecta mass for SNe Ic-BL than Taddia et al. (2019).Timescales and ejecta velocities typical of SLSNe can be reproduced over a large portion of the parameter space.However, a higher ejecta mass is required for faster spinning magnetars to keep the velocities below that of SNe Ic-BL.In contrast, a low ejecta mass is required for faster-spinning magnetars to keep the timescales below ∼ 100 days and the velocities higher than ∼ 10 000 km s −1 , providing some phenomenological justification for the mass-spin correlation found by Blanchard et al. (2020).Suzuki & Maeda (2021) find faster timescales compared to us in the L46 2D simulations, but similar timescales in their 1D simulations, and an average ejecta velocity that is roughly consistent with our model.
Figure 4 shows the peak -band absolute magnitude and peak  − colour over the same parameter grid.We find a portion of parameter space for low spin period, ejecta mass, and magnetic field which produces a transient more luminous than any previously observed superluminous supernova.However, it is unlikely such a combination (particularly low spin period and magnetic fields) can be conceived as magnetic-field amplification mechanisms such as the magnetorotational instability or Kelvin-Helmholtz instability likely amplify most typical progenitor fields to larger poloidal fields than seen in this parameter space (e.g., Reboul-Salze et al. 2021), meanwhile, the stability of magnetic-field configurations in this part of the parameter space is also questionable (e.g., Braithwaite 2009) and so magnetars in this parameter space may never materialise.We note that as we do not track gravitational-wave losses, the newly born magnetar could potentially spin-down rapidly through gravitational-wave radiation (Sarin & Lasky 2021) in this parameter space, depleting the energy reservoir to power such a luminous transient.To compound this all further, it is unknown whether stellar explosions with such small ejecta masses could harbour magnetars that are rapidly rotating but have quite weak poloidal fields in the first place.The parameter space between   = −21 and   = −23, where SLSN are, shifts to lower masses for higher spin periods, where the parameter space around   = −19, where most SNe Ic-BL are, require either a large ejecta mass (≳ 5 ⊙ ), higher spin period, or extremely high magnetic field.The parameter space where  −  < 0 mostly overlaps with the   < −21 region, showing that most SLSNe should have −0.5 <  −  < 0 at peak, which is broadly consistent with observations (Chen et al. 2023a).Suzuki & Maeda (2021) only present bolometric results, but find an increase in peak bolometric magnitude of ∼ 5 in 2D simulations and ∼ 4 in 1D simulations between the most and least energetic L46 models, while we find an increase of ∼ 4 in -band magnitude over the same parameter difference.The magnitudes are also roughly consistent with the analytical model presented in Kashiyama et al. (2016).The energy ratios from Figure 2 roughly coincide with the peak   for the parameter region where  ≲ 1, but as the spin-down time increases, a larger fraction of the emission is emitted post-peak.The  −  colours show similar behaviour to   , but show less dependence on magnetic field.We speculate that this is because the peak timescale is shorter, and thus the ejecta has had less time to cool.

Effect of Varying Braking Index
Light curve luminosity and morphology can vary significantly with variations in magnetar braking index.Figure 5 shows the bolometric luminosity, absolute -band magnitude, and absolute -band magnitude for several supernovae where only the braking index  is varied.The ejecta mass, spin-down time, and total rotational energy are fixed to 10  ⊙ , 10 6 s, and 10 52 erg, with initial magnetar luminosity calculated from Equation 4; all other parameters are the same as in Section 3.1.The timing of the light curve peak can vary by a factor of ∼ 3 and late-time luminosities by orders of magnitude, with large  peaking later and having higher luminosities at later times, although the variation in late-time luminosity asymptotes as  increases due to the exponent in Equation 3 asymptoting to −1.The peak luminosity can vary by ∼ 1 mag and is highest for  ≈ 3 in this case, although this is not true in general, and these numbers will vary depending on the energetics and diffusion time of the supernova.The requirement of constant rotational energy drives down  0 (Equation 4), leading to less energy injected at early times and more at later times, which causes the peak luminosity to decreases at large .

CASE STUDIES FOR INFERENCE
As a proof of concept, we perform inference on several different classes of supernovae to see if the model is flexible enough to recover sensible parameters for a variety of objects.First, we validate the model using a simulated SLSN.Then, we perform inference on SN 2015bn, an SLSN; SN 2007ru, a SN Ic-BL; ZTF20acigmel (better known as "the Camel"), an FBOT; and iPTF14gqr, a USSN.
Inference is performed on multiband photometry using the opensource software package Redback (Sarin et al. 2023) with the dynesty sampler (Speagle 2020) implemented in Bilby (Ashton et al. 2019;Romero-Shaw et al. 2020).We sample with a Gaussian likelihood and an additional white noise term, and sample in flux density rather than magnitude.We use the default priors in all cases (shown in Table 1) except for the explosion energy of iPTF14gqr, where the lower limit is reduced to 5 × 10 48 erg to capture the lower expected explosion energies of USSNe (Suwa et al. 2015), and the plateau temperature of ZTF20acigmel, where the observed photospheric temperature stayed at ∼ 20 000 K for the entire time it was detectable in optical/UV (Perley et al. 2021).We also sample the unknown explosion time with a uniform prior of up to 100 days before the first observation and an extinction term   with a uniform prior between 0 and 2, and use a constraint that the total rotational energy of the magnetar  rot ≲ 10 53 erg.
The fitted light curve and corner plot for the simulated SN are shown in Figure 6, the fitted light curves for the other SNe are shown in Figure 7, the input parameters for the simulated SN and recovered parameters for each SN are shown in Table 2, and the corner plots for each SNe are shown in Appendix A.

Validation on a Simulated Supernova
To test if the model could recover parameters correctly, we simulated an SLSN with  0 = 10 46 erg s −1 ,  SD = 10 6 s,  = 3,  ej = 10  ⊙ , and  = 0.1 cm 2 g −1 , with the other parameters the same as in Section 3. The supernova was placed at redshift  = 0.1 and data was generated using the Redback simulation workflow as observed by ZTF in ,  and  band for the first 200 days post explosion.
The light curve (Figure 6 (left)) is fit well by the model throughout its evolution.In the one-dimensional posteriors, only  0 is recovered to within 1, while  SD , , and  ej , are recovered to just outside 1 and  is not recovered well at all.In the two-dimensional posteriors, every parameter is recovered to within 2 except for the  - ej posterior and anything involving .The correlations between  0 ,  SD , and  also suggests that the rotational energy of the magnetar is well constrained to ∼ 10 52 erg, i.e., the injected value.The correlation between  ej and  is the main reason they are outside the 1 error region in the one-dimensional posteriors. also shows   an anti-correlation with   (see Figure A1), showing that the effect of lowering the braking index can be mimicked in some cases by incomplete gamma-ray thermalization, although this degeneracy can likely be broken by acquiring data at later times.This result shows the importance of understanding the various opacities of different types of supernovae as well as the necessity of reporting inferred opacity values.

SN 2015bn
SN 2015bn is an SLSN-I at  = 0.1136 that was first discovered by the Catalina Sky Survey on 2014 December 23.It peaked at 79 restframe days post discovery, which made it one of the slowest evolving SLSNe at the time, and had peak magnitudes of   = −22.0± 0.08 mag (AB) and   = −23.07± 0.09 mag (Vega), making it one of the most luminous as well (Nicholl et al. 2016a).Since then, it has been followed-up extensively in optical/UV/NIR, with photometry and spectroscopy (Nicholl et al. 2016a(Nicholl et al. ,b, 2018)), and polarimetry (Inserra et al. 2016;Leloudas et al. 2017), and well as in radio (Nicholl et al. 2018;Eftekhari et al. 2021;Murase et al. 2021) and X-rays (Inserra et al. 2017;Bhirombhakdi et al. 2018).SN 2015bn shows strong undulations in the light curve on a timescale of 30-50 days (Nicholl et al. 2016a;Inserra et al. 2017), and was detectable in optical/UV for more than 1000 days (Nicholl et al. 2018), although the supernova has yet to be detected in either radio or x-rays.
We import the observational data (Nicholl et al. 2016a,b) from the Open Supernova Catalog (Guillochon et al. 2017).Our model is able to fit the supernova peak very well in all bands but underestimates the IR bands in the post-peak photospheric phase and overestimates several bands after 300 days.The inferred rotational energy of the magnetar is ∼ 6 × 10 52 erg, which is close to the maximum rotational energy that can be extracted from a newborn magnetar.The gammaray opacity we find is around 10 −2.4 cm 2 g −1 , which is consistent with Vurm & Metzger (2021) at around 200-300 days, when leakage starts to become noticeable.The magnetar also shows a large spread in inferred braking indices, but excludes  = 3 at ≳ 95% confidence, meaning the vacuum dipole is not a good approximation for this object.Using Equations 1 and 4 to infer the spin period, assuming a 1.4  ⊙ neutron star with the same equation of state as Nicholl et al. (2017b), gives  0 ≈ 0.7 ms, extremely close to the mass-shedding limit.
Comparing to the results of Nicholl et al. (2017b) shows strong discrepancy between most inferred parameters.Using Equations 4, 21, and 22 with parameters from Nicholl et al. (2017b) show that they find a spin-down timescale of ∼ 105 days, much lower than our ∼ 400 days, and a magnetar rotational energy of only ∼ 10 52 erg, a factor of ∼ 5 lower than recovered by our model.This is because the magnetar needs to supply both the radiated and kinetic energy of the supernova in our model, while in the Nicholl et al. (2017b) model the ejecta velocity is input separately.We also recover a smaller ejecta mass of ∼ 6  ⊙ , compared to their ∼ 12  ⊙ , and a smaller opacity of ∼ 0.05 cm 2 g −1 , compared to their ∼ 0.19 cm 2 g −1 .
We import the observational data (Sahu et al. 2009) from the Open Supernova Catalog (Guillochon et al. 2017).The model fits most of the data well in the optical and NIR.The magnetar energy is ≲ 10 50 erg for this supernova, lower than the explosion energy, which dominates the dynamics here.This energy is also much lower than that estimated by Wang et al. (2016), and the spin-down time is a factor of ∼ 10 larger.The ejecta mass we estimate is similar to the previous nickel-powered model, but lower than the previous magnetar-powered model.The opacity we recover is ∼ 0.06 cm 2 g −1 , lower than the fixed value of 0.1 cm 2 g −1 in Wang et al. (2016).Finally, the braking index  is much higher than 3, and we can again reject vacuum dipole spin-down at ≳ 95% confidence.Our posterior on the braking index is also consistent with the magnetar spin-down being dominated by gravitational-wave radiation (Sarin et al. 2018(Sarin et al. , 2020a)).However we caution against strong conclusions based on the measurement of  due to our simplified treatment of the neutron star spin evolution.
We used the publicly available photometric data from Perley et al. (2021) for our fit.The model fits the data in all filters throughout the evolution of the object.The total rotational energy of the magnetar is around 10 51 erg, comparable to the explosion energy.Using Equations 1 and 4 to infer the spin period, again assuming a 1.4  ⊙ neutron star with the same equation of state as Nicholl et al. (2017b), gives  0 ≈ 5 ms.This value, as well as the ejecta mass, are consistent with values found for the FBOT distribution (Liu et al. 2022).The braking index has a wide distribution with two peaks, and is anti-correlated with the gamma-ray opacity, like with the simulated supernova.If the spin-down mechanism is vacuum dipole, then the gamma-ray opacity is likely ≳ 1 cm 2 g −1 , which can happen at late times for PWNe with magnetization similar to Galactic PWNe (Tanaka & Takahara 2010, 2013;Vurm & Metzger 2021).If the braking index is larger, then the gamma-ray opacity is smaller, and the magnetization must be lower, similar to SN 2015bn and SN 2017egm (Vurm & Metzger 2021).This can potentially tested by observing non-thermal emission in the late phase.

iPTF14gqr
iPTF14gqr is a USSN at  = 0.063 that was first discovered by the intermediate Palomar Transient Factory (iPTF) (Law et al. 2009) on 2014 October 14.The supernova featured a bright first peak that faded within a day, followed by a more extended light curve the rises after about ∼ 4 days.The first peak can be explained by shock cooling (De et al. 2018) and the second peak by either radioactivity (De et al. 2018) or a new born magnetar (Sawada et al. 2022), although both peaks can also be explained by interaction alone (Khatami & Kasen 2023).
We used the publicly available photometric data from De et al. (2018), although we exclude all data points within one day postexplosion, since we only claim that the second peak could be magnetar powered.The model fits the data in all filters throughout the evolution of the object.The explosion energy, which has a wider prior to account for the low explosion energy of USSNe, is constrained to ∼ 6 × 10 49 erg, while the total magnetar rotational energy is ∼ 2 × 10 49 erg.The braking index is not very well constrained; although it is consistent with vacuum dipole to within error.The median spin-down time of ≈ 7 days is slightly lower than that found by Sawada et al. (2022), although our magnetar energy is a factor of ∼ 5 higher; this is likely because our explosion energy is smaller.The ejecta mass we derive is also consistent with estimates by both De et al. (2018) and Sawada et al. (2022).

DISCUSSION AND SUMMARY
As shown by both the exploration of the parameter space (Section 3) and the case studies (Section 4), the model presented here is incredibly versatile.Three of the four case studies had previous been fit by a magnetar model (Nicholl et al. 2017b;Wang et al. 2016;Sawada et al. 2022), but each of these studies used a different model with different assumptions to model a particular supernova or class of supernovae.A versatile model allows comparisons of different populations to be done self-consistently and determine what model variations manifest in vastly different types of supernovae, as well as probe whether a continuum between these sources could exist or whether there are multiple distinct classes.
Much of the flexibility of our model comes from self-consistent dynamical evolution of the ejecta, while the addition of magnetar braking index as a parameter allows for some possible insight into the spin-down mechanism of the newborn millisecond magnetar.While inferring the spin-down mechanism for any particular supernova is difficult due to the number of possible mechanisms that could be in effect, we can constrain this parameter at a population level and see if different classes have evidence for different spin-down mechanisms.Figure 8 (top) shows the posterior probability distribution of braking index for the four case study supernova.As mentioned above,  = 3 is rejected for SN 2015bn and SN 2007ru at > 95% confidence; SN 2007ru shows a posterior that is roughly Gaussian while SN 2015bn shows a posterior that is roughly uniform at  ≳ 4. ZTF20acigmel and iPTF14gqr both show non-Gaussianity in their posteriors as well; the Camel has a double peaked distribution due to the degeneracy with gamma-ray opacity, while iPTF14gqr has a tail at high  but peaks very close to  = 3.While making definitive statements about spin-down mechanisms will require a much larger sample, this small sample already shows diversity in their inferred braking index, highlighting a potential interesting question for future studies.
All the objects studied in Section 4, including the simulated supernova, show a strong negative correlation between  0 and  SD .This shows that the magnetar rotational energy  rot , and thus the total energy budget of the supernova, can be constrained for these objects to within an order of magnitude (see Figure 8 (bottom)).The ejecta masses were all found to be similar to previous studies, except for SN 2015bn, and the spin-down times were found to differ for all models, usually by a factor of ∼ 3-10.The magnitude of these discrepancies seems to vary depending on the supernova, and a large-scale sample study on supernovae that have been previously characterized by a magnetar model (e.g.Chen et al. 2023b) is necessary to characterize the systematic difference between our model and previous models.
Each magnetar-driven transient model (e.g Yu et al. 2013;Wang et al. 2016;Kashiyama et al. 2016;Nicholl et al. 2017b;Sarin et al. 2022;Sawada et al. 2022) has slight differences in assumptions.All of these models assume vacuum dipole spin-down, although Kashiyama et al. (2016)  The kilonova models (Yu et al. 2013;Sarin et al. 2022) have different radioactive heating rates due to the r-process material in the ejecta, while the other supernova models (Wang et al. 2016;Kashiyama et al. 2016;Nicholl et al. 2017b;Sawada et al. 2022) are non-relativistic and can not be used for transients with exceptionally low ejecta masses and powerful magnetar engines.Finally, only Nicholl et al. (2017b) and Sarin et al. (2022) were originally implemented into a publicly-available Bayesian inference code (Guillochon et al. 2017;Sarin et al. 2023), although the Yu et al. (2013) model has also been implemented in Redback.
This model has a few caveats which may prevent it from properly describing certain transients.The first is that it is a one-zone, onedimensional model.This makes the treatment of the photospheric radius very simplified compared to real supernovae.Engine-driven supernovae also show hydrodynamic instabilities in multidimensional simulations (e.g.Chen et al. 2016Chen et al. , 2020;;Suzuki & Maeda 2017, 2019, 2021) which can shred the inner ejecta, causing a decrease in the effective optical depth of the ejecta, as well as affecting the timescale for non-thermal leakage.While our model can qualitatively produce the same behaviour as radiation hydrodynamic simulations, multi-dimensional effects cause the supernova peak to be earlier and more luminous.If the spin-down timescale of the magnetar is smaller than the Kelvin-Helmholtz timescale of the magnetar for neutrino emission  KH, ≲ 100 s, then baryon loading on the magnetized wind via the neutrino-driven wind can be relevant (Thompson et al. 2004) and the magnetized wind can be collimated by anisotropic and hoop stress (Bucciantini et al. 2007(Bucciantini et al. , 2008) ) and form a jet (Kashiyama et al. 2016).This model also does not self-consistently track gravitationalwave emission and how it depletes the overall rotational energy reservoir.The model also has no way to explain the bumps and undulations that have recently found in a large number of SLSNe (Hosseinzadeh et al. 2022;Chen et al. 2023b), which have been explained by both circumstellar material (e.g.West et al. 2023;Chugai & Utrobin 2023) and magnetars (Chugai & Utrobin 2022;Moriya et al. 2022;Dong et al. 2023).The model we present can only explain a smooth light curve with minimal fine structure; however, this structure is likely connected to small mass ejections or binary interactions in the final few years of the life of the progenitor and may not be strongly connected to the power source of the supernova.Furthermore, The SED used in our model has a simplified treatment of line blanketing that is calibrated to SLSNe (Nicholl et al. 2017b), and may not be accurate for other transients.Although it is possible to switch the SED to a blackbody if line blanketing is not expected to be strong, the SED is also not a good approximation deep into the nebular phase on the supernova, since the emission will start to be dominated by line emission instead of photospheric emission (e.g., Schulze et al. 2023).
Although the magnetar-driven supernova model is versatile enough to fit light curves of many different supernovae, the best way to determine whether a magnetar is really the power source is to compare the nebular spectra, polarization, and non-thermal emission from the supernovae with different models.Within the magnetar model, this emission is form the PWN or its interaction with the supernova ejecta.If the ejecta is optically thin, the synchrotron and inverse Compton emission from the PWN can leak through and be detected (Kotera et al. 2013;Metzger et al. 2014;Murase et al. 2015;Omand et al. 2018), while if the ejecta is optically thick, these photons will be absorbed and change the temperature and electronic state of the ejecta, both giving detectable signals (Chevalier & Fransson 1992;Omand et al. 2019;Omand & Jerkstrand 2023).However, interaction with circumstellar material can also produce non-thermal emission, polarization, and spectra with high ionization lines, therefore detailed modeling is necessary to make any strong conclusions.
In this work, we present a more flexible, inference-capable, publicly available model for the light curves of magnetar-driven supernovae.This model can be applied to any transient where the dominant power sources are either a magnetar or 56 Ni, as long as that energy release is roughly spherically symmetric.The main changes from previous models are the coupling of the magnetar energy injection to the kinetic energy of the supernova and the addition of the magnetar braking index as a free parameter, allowing the exploration of non-vacuum-dipole spin down.We show that the model can reproduce the basic properties on several phenomenologically different supernovae, and also fit four different types of supernovae, retrieving parameters consistent with works using separate models.This model will allow us to explore the full diversity of these supernovae, better characterize these supernovae from just their light curves, and make better predictions future multiwavelength emission to better test the magnetar-driven scenario.

Figure 1 .
Figure 1.Initial pulsar luminosities  0 (top) and spin-down times  SD (bottom) for different initial pulsar rotation periods and magnetic fields.

Figure 2 .Figure 3 .
Figure 2. Ratio of kinetic to radiated energy for supernovae with varying ejecta mass and  0 = 1 ms (left) and  0 = 3 ms (middle) and with varying spin period and  ej = 10 ⊙ (right).The black lines indicate contours of constant  kin / rad (1,10, and 100), while the magenta lines indicate contours of constant  =  SD / dif (0.1, 1, and 10).The black symbols indicate the parameters used in Suzuki & Maeda (2021), with the circles representing the L46 models and the squares representing the L48 models.

Figure 5 .
Figure 5. Bolometric luminosity (left), absolute -band magnitude (middle), and absolute -band magnitude (right) for several supernovae where only the braking index  is varied.

Figure 6 .
Figure 6.Fitted light curve (left) and posteriors of key parameters (right) for the simulated SLSN.The solid lines in the light curve plot indicate the light curve from the model with the highest likelihood, while the shaded area indicates the 90% credible interval.The orange dots and lines in the posterior indicate the injected parameters.

Figure 7 .
Figure 7. Fitted light curves for the four supernovae from our case studies.The solid lines indicate the light curve from the model with the highest likelihood, while the shaded area indicates the 90% credible interval.

Figure 8 .
Figure 8. Posterior distributions of the magnetar braking index (top) and rotational energy (bottom) for the four supernovae from our case studies.The median values are indicated by the blue vertical lines within the distribution.The black dashed line indicates vacuum dipole spin-down ( = 3).
and Sarin et al. (2022) also add a spin-down component for gravitational waves from magnetic deformation.For gamma-ray thermalization, Yu et al. (2013) uses a constant gammaray thermalization efficiency, Kashiyama et al. (2016) models the non-thermal emission and deposition using several additional free parameters, while Wang et al. (2016), Nicholl et al. (2017b), and Sarin et al. (2022) use the same prescription we do.Nicholl et al. (2017b) and Sawada et al. (2022) do not treat the effect of the magnetar on the ejecta dynamics, while Wang et al. (2016) uses a simplified treatment.

Figure A1 .
Figure A1.Corner plot for the simulated SLSN.Injected values are shown in orange.

Table 2 .
Injected parameters for the simulated supernovae and median inferred parameter values and 1 uncertainties for the simulated supernova and the four supernovae from the case studies.