Magnetar Flare-Driven Bumpy Declining Light Curves in Hydrogen-poor Superluminous Supernovae

Recent observations indicate that hydrogen-poor superluminous supernovae often display bumpy declining light curves. However, the cause of these undulations remains unclear. In this paper, we have improved the magnetar model, which includes flare activities. We present a systematic analysis of a well-observed SLSNe-I sample with bumpy light curves in the late-phase. These SLSNe-I were identified from multiple transient surveys, such as the Pan-STARRS1 Medium Deep Survey (PS1 MDS) and the Zwicky Transient Facility (ZTF). Our study provides a set of magnetar-powered model light curve fits for five SLSNe-I, which accurately reproduce observed light curves using reasonable physical parameters. By extracting essential characteristics of both explosions and central engines, these fits provide valuable insights into investigating their potential association with gamma ray burst engines. We found that the SLSN flares tend to be the dim and long extension of the GRB flares in the peak luminosity versus peak time plane. Conducting large-scale, high cadence surveys in the near future could enhance our comprehension of both SLSN undulation properties and their potential relationship with GRBs by modeling their light curve characteristics.


INTRODUCTION
Superluminous supernovae (SLSNe) have been a hot topic in supernova research for nearly two decades since the first discovered events (e.g., SN 2005ap; Quimby et al. 2007). SLSNe are characterized by their high luminosity, which can be up to tens to hundreds of times brighter than typical supernovae (see Gal-Yam 2019, for a review). SLSNe can be divided into two categories based on their spectra: hydrogenrich (SLSN-II) and hydrogen-poor (SLSN-I). It is widely believed that the high luminosity of SLSN-II is primarily attributed to the interaction between the supernova ejecta with their massive hydrogen-rich circumstellar medium (CSM) (Smith & McCray 2007;Chevalier & Irwin 2011). The total energy radiated by a typical SLSN-I is approximately a few times 10 51 ergs, if all of this energy is provided by the decay of 56 Ni, it would require at least several solar masses of 56 Ni. However, it is very difficult for a regular supernova nucleosynthesis process to produce so much 56 Ni (Umeda & Nomoto 2008). The conventional model of radioactive power is facing a significant challenge due to the exceptionally high luminosity exhibited by SLSN-I 1 . There are currently two main scenarios for the energy source of SLSN-I, namely the CSM interaction scenario (Chevalier & Irwin 2011;Ginzburg & Balberg 2012;Chatzopoulos et al. 2012) and the magnetar-powered scenario (Kasen & Bildsten 2010;Woosley 2010). The former scenario is simi-1 LCs of some slowly evolving SLSNe could be explained by radioactive decay of massive amounts of 56 Ni synthesized in a pair-instability supernova explosion (e.g., Gal-Yam et al. 2009;Gutiérrez et al. 2022). However, no SLSN has been convincingly shown to be consistent with a radioactivepowered model. Even those that have a similar decay rate as 56 Co decay are too blue and rise too quickly. The vast majority of rapidly evolving SLSNe-I are unlikely to be primarily powered by radioactive decay. lar to SLSN-II, except that the SN is surrounded by a massive hydrogen-poor CSM. The latter scenario considers energy injection from the spin-down of a newborn rapidly rotating magnetar. Numerous prior studies have demonstrated that these two scenarios are capable of explaining the general shape of the photometric light curves (LCs) of SLSNe-I (e.g., Chatzopoulos et al. 2013;Inserra et al. 2013;Wang et al. 2015;Prajs et al. 2017;Liu et al. 2017;Nicholl et al. 2017;Chen et al. 2023), which makes it difficult to distinguish one scenario from the other. Further photometric observations are necessary to discriminate between the energy source models. Currently, time-domain surveys are gathering substantial amounts of well-observed SLSNe-I. The magnetar-powered model or CSM interaction alone cannot explain undulations and occasional significant post-peak bumps observed in some SLSNe-I LCs (e.g., Nicholl et al. 2016a;Inserra et al. 2017;Yan et al. 2020;Gomez et al. 2021;Fiore et al. 2021;West et al. 2023;Gutiérrez et al. 2022;Tinyanont et al. 2022;Zhu et al. 2023). Recently, systematic studies of the prevalence and properties of post-peak bumps of SLSNe have found that these fluctuation structures may be relatively common (Hosseinzadeh et al. 2022;Chen et al. 2023).
The identification of post-peak bumps in certain SLSNe-I may provide novel understanding into the energy source mechanisms of these SNe. Liu et al. (2018) has shown that a multiple ejecta-CSM interaction model can explain the undulating bolometric LC of iPTF 15esb and iPTF 13dcc. The absence of relatively narrow emission lines in most spectral readings of SLSNe-I is often seen as a major challenge to the CSM-interaction as the primary power source.
In the magnetar-powered model, the LCs of SLSNe are determined by the behavior of the central engine over time. If the spin-down of the magnetar is mainly through magnetic dipole radiation, then a relatively smooth light curve can be expected. In fact, magnetars have intermittent and violent energy burst activities, which is an important observational feature of magnetars in the Milky Way (see Kaspi & Beloborodov 2017, for a review). Furthermore, in about half of Swift gamma-ray bursts (GRBs) X-ray afterglows, a large number of rapidly rising and falling X-ray flares have been discovered (Burrows et al. 2005;Falcone et al. 2007;Wang & Dai 2013;Yi et al. 2016). This strongly suggests that the GRB engines have undergone numerous intermittent and energetic activities. It is proposed that a millisecond spindown magnetar could simultaneously power long GRBs and SLSNe (e.g., Margalit et al. 2018). They both have a similar host preference in metal-poor dwarf galaxies with relatively high star formation rate (e.g., Lunnan et al. 2014;Chen et al. 2015;Perley et al. 2016;Chen et al. 2017;Schulze et al. 2018). Besides, the spectra of some SLSNe-I in their nebular phase are remarkably similar to those of GRB-SNe (e.g., Nicholl et al. 2016b;Jerkstrand et al. 2017). Given the significant correlation between GRBs and SLSNe (as note by Metzger et al. 2018;Liu et al. 2022), it is reasonable to assume that flare activity may also be caused by the central magnetar of SLSNe. These intermittent flare activities may reheat the SNe ejecta in later stages, and then show undulations characteristics on the LC .
This paper focuses on the origin of late-time undulations and whether the magnetar model with flare components can replicate the observed bumpy declining light curves (LCs) of SLSNe-I. To achieve this, we have updated the magnetar model by incorporating the Markov chain Monte Carlo (MCMC) technique and conducted a systematic analysis using well-sampled late-phase LCs of SLSNe-I, including multi-color LC fits and parameter comparisons. The paper is structured as follows: Section 2 describes a sample of wellobserved SLSNe-I with post-peak bumps; Section 3 provides a brief description of our model and fitting process for the LCs; Section 4 presents statistics on fitting parameters and compares them to GRB flares; finally, in Section 5, we draw conclusions from our findings and discuss their implications.

SAMPLE SELECTION
Previously, Hosseinzadeh et al. (2022) found that the majority (44-76%) of their SLSNe-I sample possess late-time LC fluctuations. Similar undulations were detected in 34%-62% SLSNe-I discovered in ZTF Phase-I Survey (Chen et al. 2023). However, these work mostly focused on observational statistics, whereas the existing theoretical study of LC fluctuations mostly focused on individual objects (e.g., Fiore et al. 2021;Chugai & Utrobin 2022;Tinyanont et al. 2022;West et al. 2023). Here, we perform systematic fits on SLSNe-I LCs with bumpy features under the magnetarpowered model. To investigate the undulations of SLSNe-I, we need adequate observational data on their post-peak bumps. We collected publicly available data on SLSNe-I in the Open Astronomy Catalog (OAC) 2 and the Bright Transient Survey explorer (BTS, Fremling et al. 2020;Perley et al. 2020) catalog 3 . Our sample includes all SLSNe-I that meet the criteria listed in Table 1. Specifically, We performed the selection with the following three steps: 1. This object has been spectroscopically classified as a SLSN-I.
2. The light curve of SLSN must be well observed with a distinct primary peak in at least two filters, accompanied by at least five data points for both the rising and declining phases.
3. The SLSN late-time LCs deviate from the smooth decline trend predicted by the magnetar spin-down model and feature in one or more additional post-peak "bumps", which contain a complete rising and declining structure (i.e., an obvious second or third peak) and should be identified with residuals larger than 3σ. Table 1. Criteria for selecting SLSNe-I with late-time bumps.
Step Critiria Sample Size 1 Spectroscopically identified as SLSNe-I 182 2 Main peak well sampled & 80 no large LC gaps 3 >3σ deviation from the magnetar model & 5 complete second/third peak(s) In addition, SLSNe-I with large gaps or deficient rising phase in their LCs (e.g., SN 2017gci) are excluded from our sample. Following the criteria, we finally selected 5 SLSNe-I 4 . Their basic information and references are listed in Table 2. For the three ZTF objects of our sample, we obtained their apparent magnitude using the python-based package HAFFET (Yang & Sollerman 2023), which can request data from the ZTF Forced Photometry Services (Masci et al. 2019) directly. In our sample, three SLSNe only contain observations in two bands. To avoid introducing additional systematic errors, we fit their multi-color LCs instead of a bolometric one, which is also able to unveil the color information. Our model assumes that the continuous injection of energy determines the underlying profile of SLSN LCs. The undulations in the LC are caused by flare activities from the central magnetar. Initially, we examine when and how SLSN LCs deviate from the smooth magnetar model. For each SN in our sample, we subtract the best-fitting model from the observed LC to identify characteristics of late-time excess. Subsequently, we will establish a connection between these post-peak bumps and magnetar flare activities.

The underlying profile of the SLSNe light curve
In the previous magnetar-powered model, the typical assumption is that the spin-down luminosity of the magnetar follows the magnetic dipole radiation as (e.g., Chatzopoulos et al. 2012;Inserra et al. 2013;Nicholl et al. 2017) where L sd,i ≃ 10 47 erg s −1 P −4 i,−3 B 2 p,14 is the initial spindown luminosity and t sd ≃ 2 × 10 5 sP 2 i,−3 B −2 p,14 is the characteristic spin-down timescale. P i and B p are the initial spin period and polar magnetic strength of the magnetar, respectively. Here we take the conventional notation Q x = Q/10 x in cgs units. This spin-down luminosity of the magnetar inflates a nebula of relativistic particles and radiation inside the cavity evacuated by the expanding supernova ejecta.
Since the energy from the magnetar may be in the form of high-energy photons, we can incorporate the gamma ray leakage effect into the standard supernova diffusion equation to calculate the bolometric luminosity as (Arnett 1982) where the diffusion time t diff = (2κM ej /βcv ej ) 1 2 , determined by ejecta mass M ej , velocity v ej , and optical opacity κ. β is a constant integration ∼ 13.8 and c is the speed of light. The factor 1 − e −At −2 accounts for high-energy photons, and A = 3κ γ M ej /4πv 2 ej (Wang et al. 2015), κ γ is the opacity of gamma rays.
To calculate the multi-band LC of a SLSN, we assume a photospheric temperature as ) where σ SB is the Stefan -Boltzmann constant and T floor is the photospheric temperature. The green, red, and orange solid lines represent the best baseline fitting results in g, r, and i band, respectively. 2σ interval is shown as translucent color range. Residuals with respect to the baseline are shown in the lower panels. Data inside the gray shaded region were used to determine the baseline. The relatively larger error range of PS1-12cil in g and r band is attributed to the large sampling interval.
In order to obtain the baseline parameters properly, we only utilized observations of the main peak in their LCs to fit the baseline (marked as the shaded grey region in Figure 1). The model fittings were performed using the MCMC technique based on the python-based package emcee (Foreman-Mackey et al. 2013). The typical number of iterations for the MCMC algorithm is 5000. All the free parameters used in this model are listed in Table 3. The maximum likelihood approach was adopted to obtain the best fitting parameters, which are listed in Table 4. We present the best-fitting baseline with residuals for each SLSN in Figure 1, where the underlying LC profile of each SLSN is well described by the magnetar model. With respect to the basic model, the late-time deviations are rather distinct for all 5 sources, which clearly requires an additional physical process. Moreover, the luminosity excess of each source appears almost simultaneously in different filters, suggesting that the bump is wavelength-independent.

Light-curve undulations and magnetar flare activities
Besides magnetic dipole radiation, another way to release the rotation of a magnetar is through erupting magnetic bubbles due to the wind-up of the seed field (Kluźniak & Ruderman 1998;Ruderman et al. 2000). This process naturally creates an unpredictable central engine that is significant in understanding the prompt emission of GRBs and X-ray flares (Dai et al. 2006). In nearly half of all GRBs, bright X-ray flares are detected (Burrows et al. 2005;Falcone et al. 2007;Wang & Dai 2013;Yi et al. 2016). This strongly suggests that the GRB central magnetar have undergone numerous intermittent and energetic activities. By considering the strong correlation between GRBs and SLSNe Metzger et al. 2018;Liu et al. 2022), it is reasonable to assume that the central magnetar of SLSNe could also trigger flare activities, leading to undulations in their LCs due to delayed release of flaring energy.
We describe the energy release rate resulting from the central magnetar flare using a formula commonly used to fit GRB X-ray flare LCs (Yi et al. 2016) where α 1 , α 2 and ω are the structure parameters reflecting the sharpness and smoothness of the flares, which have relatively small impact on the final LC shapes . In this work, we adopted general fixed values with α 1 = −2, α 2 = 2.5 and ω = 5. L flare,pk and t pk are free parameters representing the flare energy injection and peak time. However, unlike in the case of GRB flares, the energy and time of flare activities cannot be straightforwardly compared to the observed post-peak bumps of SLSNe-I. The main reason is that the flare energy is likely trapped in the optically thick SN ejecta and re-radiated by diffusion. The energy released from the magnetar flare activities can reheat SN ejecta and then cause the "bump" feature appearing in the LCs. By substituting the expression of L flare (t) + L sd (t) into Eq.(2) to replace L sd (t), we can obtain the LC after being affected by the flare activities of central magnetar.
By adding additional flare components to the best-fitting baseline, we successfully reproduced the multi-color LCs of the 5 sources. The final fitting results are presented in Figure 3, where two flares are introduced to explain the LCs of SN 2019lsq and one for other objects. The corresponding parameters are listed in Table 5. For clearer demonstration, we also plot the bolometric LCs of two sources (i.e., SN 2019stc and PS1-12cil) in Figure 4, along with the underlying energy injection and the flare component respectively. The bolometric LCs were constructed using SuperBol 5 (Nicholl 2018) based on the observed fluxes in available filters. We simply fit their spectral energy distribution (SED) assuming blackbody radiation 6 . Time dilation and the Galactic extinction were corrected. We succeed in reproducing the bumpy bolometric LCs by performing the same fitting procedure adopted for multi-color LCs. The best parameters are consistent with those from the multi-color fitting. NOTE-The last three parameters are highly dependent on the light curve features of the specific object, where t flare,rise is the theoretical estimation rather than an observable. Therefore, the chosen prior ranges are general.

RESULTS
In this section, we present the quantitative analysis of the fitting parameters obtained for the basic magnetar and flares, and also compare the flare properties with those of GRBs. The three key parameters for the underlying magnetar are M ej , P i and B p . In our sample, the ejecta mass ranges in (2.22 − 18.03) M ⊙ , the initial period are distributed in (1.73 − 6.24) ms, and the magnetic field scatter between (0.25 − 3.09) × 10 14 G. We plot the key parameters of our sample in two parameter planes (the M ej and B p versus P i ) in Figure 2, where 31 SLSNe-I are also shown for comparison. The baseline parameters of our sample roughly lie within the typical range of SLSNe-I, indicating that the basic energy injection of the undulated sources in our work has no specificity.
In terms of flares, we extracted six flares from the 5 selected SLSNe-I. All the flares were injected (3.4 − 7.4) ×  10 6 s after the explosion. The rise time distributes in a range of (1.8−7.3)×10 6 s, and mainly within (1.8−3.5)×10 6 s. Their peak luminosity ranges in (0.85 − 43) × 10 43 erg s −1 . NOTE-t flare,rise represents the rise time of a flare luminosity from injection to the peak, which can be converted to the flare peak time (t pk ) by adding the t0.
In general, the SLSN flares in our sample last much longer and are much dimmer than the GRB flares (with the peak times and duration mainly between 100 − 1000 s (Yi et al. 2016)). We plot the flare parameters (L flare,pk and t pk 7 ) of our sample with 200 GRB flares from Yi et al. (2016) in Figure 5. They found that the GRB flare parameters follows an empirical relation of log L flare,pk = 51.73 − 1.27 log t pk .
The study conducted by  reveals that the three flares of SN 2015bn are situated within the 2σ range of the relation established in the context of GRB flares (represented by blue stars in Figure 5). A similar result was found in our work. After adding flares of our sample, we refit the GRB and SLSN flare parameters and found that the relation still stands, which can be expressed as log L flare,pk = 51.54 − 1.18 log t pk .
This relation has a stronger Pearson coefficient reaching −0.85, suggesting that the more luminous the flare, the later it appears. However, we could not find such a relation in the  SLSN sample only, which may be due to the small sample size.

CONCLUSIONS AND DISCUSSION
We have presented the systematic analysis of a sample of SLSNe-I with late-time undulations under the central magnetar assumption. We employed a semi-analytical magnetar model with additional flare activities to fit their multi-color LCs, as well as two available bolometric LCs, and compared the flare properties with GRB flares. Our main conclusions are as follows; 1. The undulating LCs of the 5 selected SLSNe-I can be explained by additional flare activities from a central spin-down magnetar. On the one hand, the underly-ing LC profile of these objects is consistent with the magnetar-driven baseline, and the corresponding parameters (P i , M ej and B p ) fall within the typical range observed for SLSNe-I (see Figure 2), suggesting no peculiar features existing in their basic energy sources. On the other hand, the distinct late-time bumps in their multi-color LCs can be accurately replicated by invoking magnetar flares.
2. Six flares were extracted from the LCs of the 5 selected SLSNe-I. The injection of flares started between day 39 and day 86 after the explosion, and it took 21 to 84 days to reach the peak luminosity, which ranged in (0.85 − 43) × 10 43 erg s −1 . Flare parameters of SN 2015bn (blue star) are taken from . The solid black line represents the correlation between the peak luminosities and peak times taking both the GRB and SLSN flares into account. 3σ region of the fitting line is shown as the grey shadow. The small box on the right zooms in the SLSNe flares. The Pearson correlation coefficient "r" is shown.
3. SLSN flares tend to be systematically dimmer and longer than GRB flares. The strong correlation between peak luminosity and peak time in GRB flares remains evident and is even stronger (r = −0.85) when including the SLSNe-I sample. This is reasonable given the intrinsically weaker magnetic field in SLSNe-I compared to GRBs. However, we do not find a similar correlation in the SLSN sample alone, possibly due to the small sample size.
In previous literature, different models have been discussed to explain the late-time undulations. Within the framework of magnetar models, Chugai & Utrobin (2022) proposed that bumps are caused by the magnetar dipole field enhancement several months after the explosion. Moriya et al. (2022) suggested that the light-curve bumps are caused by variations of the thermal energy injection from magnetar spin down. However, Chugai & Utrobin (2022) only provide an explanation for a single bump and Moriya et al. (2022) have predicted an increase in photospheric temper-ature which was not observed in SN 2020qlb. In addition to the magnetar model, other models also include the interaction model (e.g., Inserra et al. 2017;Fiore et al. 2021;Li et al. 2020;West et al. 2023), the sudden drop of the ejecta opacity (Metzger et al. 2014), neutron stars colliding with binary companions (Hirai & Podsiadlowski 2022). Undulations may also be caused by an asymmetry in the geometry of the ejected material (Kaplan & Soker 2020). In this paper, we believe that the late-stage flare activity of the central magnetar can well explain the undulation in the LC.
In general, the magnetar assumption accommodates the systematic diversity of the massive progenitors, allowing for a wide range of magnetic fields to exist in magnetars with similar behaviors. This aligns with the correlation between late-time undulations of SLSNe-I and GRB flares in the L flare,pk verses t pk plane, providing supporting evidence that they share a common power origin. The linearity between peak time and luminosity over many orders of magnitude demonstrates that a dimmer flare is likely to peak at a later time, with a consistent radiated energy of ∼ 10 51 erg. This suggests that there may be regular pulsive activities of these newborn magnetars. If this is the case, the flares of SLSNe-I may also be produced through magnetic reconnection processes similar to solar flares (Shibata & Magara 2011), resulting in shock heating, particle acceleration, and energy ejection.
Although we have drawn conclusions, it is important to note that our findings are based on a limited sample size. To gain further insight into the power source of SLSNe and investigate any potential relation between SLSN and GRB, more observations of well-sampled undulated phases are necessary. Fortunately, upcoming telescopes and surveys such as Legacy Survey of Space and Time (LSST), ZTF phase II, and James Webb Space Telescope (JWST) will likely discover numerous SLSNe with clearly defined light curves in the near future. This will allow for a more comprehensive understanding of their explosion mechanisms and corresponding progenitors.
We thank Zongkai Peng and Zhihao Chen for helpful comments and discussions, and an anonymous referee for important suggestions. This work is supported by the National Natural Science Foundation of China (Projects 12021003), the National SKA Program of China (2022SKA0130100), and the National Key R&D Program of China (2021YFA0718500). The ZTF forced-photometry service was funded under the Heising-Simons Foundation grant #12540303 (PI: Graham).