Concise Spectrotemporal Studies of Magnetar SGR J1935+2154 Bursts

SGR J1935+2154 has truly been the most prolific magnetar over the last decade: it has been entering into burst active episodes once every 1–2 yr since its discovery in 2014, it emitted the first Galactic fast radio burst associated with an X-ray burst in 2020, and it has emitted hundreds of energetic short bursts. Here, we present the time-resolved spectral analysis of 51 bright bursts from SGR J1935+2154. Unlike conventional time-resolved X-ray spectroscopic studies in the literature, we follow a two-step approach to probe true spectral evolution. For each burst, we first extract spectral information from overlapping time segments, fit them with three continuum models, and employ a machine-learning-based clustering algorithm to identify time segments that provide the largest spectral variations during each burst. We then extract spectra from those nonoverlapping (clustered) time segments and fit them again with the three models: the cutoff power-law model, the sum of two blackbody functions, and the model considering the emission of a modified blackbody undergoing resonant cyclotron scattering, which is applied systematically at this scale for the first time. Our novel technique allowed us to establish the genuine spectral evolution of magnetar bursts. We discuss the implications of our results and compare their collective behavior with the average burst properties of other magnetars.


INTRODUCTION
High-energy transient activity from magnetars falls into three general categories: short and intermediate bursts, and Giant Flares (GFs).The most frequently observed events are the short bursts, with durations ranging from a few milliseconds to a few seconds, peaking at ∼ 0.1 s; their energies approach 10 41 ergs (Göǧüş et al. 2001).GFs cover the other end of the magnetar burst properties with much harder spectra and longer durations.These are the most energetic events, releasing over 10 44 ergs in several minutes (Hurley et al. 1999;Palmer et al. 2005).Apart from their longer durations, they also exhibit a unique morphology: a spectrally hard initial short spike, followed by a longer tail, that oscillates at the spin frequency of the parent neutron star.Finally, intermediate events exhibit energetics and durations in between short bursts and GFs.Their durations range from a few seconds to a few tens of seconds with energies extending up to ∼ 10 42 ergs (Ibrahim et al. 2001;Kouveliotou et al. 2001).
The run-of-the-mill magnetar bursts are generally attributed to local or global yielding of the solid neutron star crust under the influence of enormous internal and external magnetic pressure, along with excitation of modes by instabilities in the magnetosphere (Thompson & Duncan 1995, 2001).This idea was put into action by Lander (2023) who showed that the build-up of elastic stress in the crust results in its failure and in energy release.Another possible process that can generate magnetar bursts is magnetic field line reconnection (Thompson & Duncan 1995;Lyutikov 2003).
Hard X-ray spectral analyses of magnetar bursts are crucial towards a better understanding of the physical mechanisms responsible for these intriguing phenomena.Studies so far have shown that magnetar bursts can be modeled almost equally well with a thermal model as with the sum of two black bodies (BB+BB), or a non-thermal model, such as a power law with a high-energy exponential cutoff (COMPT) (van der Horst et al. 2012;Lin et al. 2012).This poses a puzzle in the understanding of the underlying mechanism of the observed bursts, as follows.According to the crustal fracturing scenario (Thompson & Duncan 1995), a fireball by the plasma of trapped e − −e + pairs and photons forms in the closed magnetic field lines by Alfven waves released from the crust following the cracking.Therefore, thermal radiation could be expected from such regions in quasi-equilibrium.On the other hand, the observed magnetar synchrotron-like non-thermal radiation spectra might be an indication of particle interactions via magnetic reconnection.Recently, Yamasaki et al. (2020) studied the spectral modification of extremely energetic magnetar flares by resonant cyclotron scattering and showed that the scattering process may alter the emerging radiation from these events.In this model, photons are emitted via a mechanism that depends only on the temperature of the fireball near the neutron star surface interact with magnetospheric particles, resulting in significant changes in the emission spectrum.
Time-resolved X-ray spectroscopy of magnetar bursts is an important probe to reveal spectral evolution throughout their highly-complex emission episodes.One of the most comprehensive timeresolved burst spectral studies was performed by Israel et al. (2008) on the SGR 1900+14 bursts observed with the Burst Alert Telescope (BAT) and X-ray Telescope (XRT) on board the Neil Gehrels Swift Observatory (hereafter Swift) in March 2006.Although the BB+BB model was at the forefront of this study including over 700 extracted spectra, a significant number of the spectra were also fitted with non-thermal models.Detailed analyses based on the BB+BB model revealed significant spectral evolution, even during short magnetar bursts.Younes et al. (2014) also performed time-resolved spectral analysis of the 63 brightest bursts from SGR J1550−5418 observed with the Gamma-ray Burst Monitor (GBM) on the Fermi Gamma-ray Space Telescope (hereafter Fermi ) and demonstrated flux dependent variations between temperature obtained with the BB+BB model and the area of the emitting region.
Time binning for the extraction of spectra in time-resolved spectroscopy of bursts, however, is quite arbitrary.In earlier studies, the subsequent spectra were usually obtained from time intervals determined based upon the signal-to-noise ratio of its light curve.In other words, the time segments are not determined by taking into account the spectral changes in advance, but the burst is divided blindly into time segments to contain a certain number of photons to ensure acceptable statistics for the spectral analysis.Under such circumstances, it would not be possible to elucidate the true spectral evolution of the observed burst emission.
One possible way to overcome the problem of arbitrariness in time-resolved spectroscopy is to employ sequential spectra extracted from the shortest possible time intervals (in terms of counts) and allow these segments to overlap subsequently.This way, it would naturally be possible to reveal real spectral evolution, which in turn could help uncover more dominant underlying mechanism.This sliding time window approach is a methodology often used in timing analysis, e.g., searching for timedependent burst oscillation behaviors (see e.g., Strohmayer et al. 1998).In this study, we apply a similar approach for the first time in the spectral analysis of SGR bursts; namely, "overlapping timeresolved spectroscopy", using the brightest bursts from SGR J1935+2154.We subsequently apply a machine learning based clustering algorithm to form non-overlapping time intervals with significant spectral variations and perform "clustering-based" time-resolved spectral analysis of these bursts.Hence, the resulting time segments are expected to precisely reveal the burst spectral evolution.
To perform spectral modeling of SGR J1935+2154 bursts, we employ commonly-used BB+BB and COMPT models.In addition, we also employ the model of a modified black body (MBB) spectrum (Lyubarsky 2002) that undergoes resonant cyclotron scattering (RCS), the combination of which is applied systematically to the SGR burst spectral analysis at this scale for the first time.The MBB-RCS model was developed by Yamasaki et al. (2020) to account for the magnetospheric effects on thermal emission in the context of magnetar flares.Lyubarsky (2002) had previously suggested extreme magnetic fields of magnetars could modify the emerging radiation from its surface, incurring significant deviations from a Planckian that result in a flat photon spectrum at low energies.However, the pure MBB model can underestimate the spectra at high energies by not taking the magnetospheric scatterings into account.The MBB-RCS model (Yamasaki et al. 2020) assumes that the photons emitted from the active burst region escape to infinity after just one resonant scattering by magnetospheric charges, thereby accounting for tails in the observed spectra at high energies.Testing the model with energetic flares of SGR 1900+14 resulted in a good agreement with the spectra of intermediate flares, but not with the giant flare observed in 1998.
SGR J1935+2154 was discovered on 2014 July 5, when a short burst triggered Swift/BAT.Followup observations with Chandra and XMM-Newton revealed its spin period (P ∼ 3.24 s) and spindown rate ( Ṗ = 1.43 × 10 −11 s s −1 ), corresponding to a magnetic field strength of B ∼ 2.2 × 10 14 G, and thus confirming the magnetar nature of the source (Israel et al. 2016).Since its discovery, SGR J1935+2154 has been the most prolific transient magnetar ever observed: it is burst-active almost annually including multiple (in the range of thousand) short bursts (Lin et al. 2020a,b), an intermediate flare on 2015 April 12 (Kozlova et al. 2016), and a burst forest on 2020 April 27 (Kaneko et al. 2021).Just hours after the burst forest, SGR J1935+2154 emitted an X-ray burst (e.g., Mereghetti et al. 2020) coincident with a Fast Radio Burst (FRB; CHIME/FRB Collaboration et al. 2020;Bochenek et al. 2020), which was the first Galactic detection of these events.This coincidence is the first indicator that magnetars residing in distant galaxies may also be the origin of FRBs (Petroff et al. 2019(Petroff et al. , 2021)).
This paper is organized as follows: In Section 2, we introduce instrument, data, the SGR J1935+2154 bursts that we studied, and their spectral data extraction process.In Section 3, we explain the steps in the spectral analysis: first overlapping and then clustering-based time-resolved spectroscopy.The results are presented and discussed in Section 4.

INSTRUMENT & SGR J1935+2154 BURST OBSERVATIONS
Fermi has been providing an enormous amount of data that allows for studying a wide range of gamma-ray transient events over the last 15 years.It carries two instruments: the Gamma-ray Burst Monitor (GBM; ∼ 8 keV−40 MeV) and the Large Area Telescope (LAT; ∼ 20 MeV−300 GeV).GBM consists of 12 sodium iodide (NaI) detectors (∼ 8 keV−1 MeV) and two bismuth germanate (BGO) detectors (∼ 200 keV−40 MeV) (Meegan et al. 2009).Since the bulk of emission from magnetar bursts is typically seen in ≲ 200 keV, we only used the data of NaI detectors for this study.We employed Continuous Time-Tagged Event (CTTE) data of GBM, which provides the finest time (2 µs) and energy (128 channels) resolutions.Our investigations were performed in the 8−200 keV energy band with 4 ms minimum time resolution.For each burst, we included data collected with the three brightest NaI detectors1 with the detector-to-source angle being less than 60 o .Additionally, we excluded detectors if they are partially or fully blocked by other parts of the spacecraft as obtained using the GBMBLOCK software provided by the GBM team.
We selected 51 SGR J1935+2154 bursts observed with Fermi -GBM between 2014 and 2022 for our time-resolved spectral investigations, based on the results of time-integrated analyses (Lin et al. 2020a,b, Lin et al. in preparation).In particular, we chose the bursts that contain at least 2400 background-subtracted counts to ensure that they have enough statistics for time-resolved spectral analysis (See Appendix A).Our investigations are done in two stages for each burst: First, we define the overlapping time segments with which we obtain spectral parameters and use a machine learning clustering algorithm to obtain "change points" for the parameter evolution.Then, we analyze spectral data extracted from the (non-overlapping) time intervals between these change points to better characterize spectral evolution.
In the first stage of our time-resolved spectral studies, we required each time segment to overlap with the previous time segment by 80%.In other words, the subsequent time segment does not start from the end of the previous one, but from the time that is one fifth of the previous segment (see Figure 1 for an example).We also required a minimum of 1200 background-subtracted counts in each time segment (See Appendix A).This way, burst spectral parameters are expected to be well-constrained and statistically reliable throughout each burst.We note that for each burst, we calculated duration2 via Bayesian Block method optimized for photon counting time series (Scargle et al. 2013), and started the first time segment from the beginning of Bayesian Block duration.
Time segments starting before the peak of a burst are prone to end at the peak although they are shifted by 20% of the time length of their previous ones due to the fact that flux rise timescale .Light curve of an SGR J1935+2154 burst detected at 2021 January 29 10:35:39.918UTC as seen with the brightest detector (n5).The vertical dashed lines show Bayesian Block duration start and end times, respectively.The red horizontal lines represent the 22 overlapping time segments with each subsequent segment having an overlap of 80%.We note that the vertical values corresponding to the red lines are arbitrary, and chosen only for display purposes.
in short magnetar bursts is usually shorter than that of decay (Göǧüş et al. 2001).This yields an accumulation of time segments at the beginning of the bursts.In those cases, we avoided the accumulation of time segments by reducing the overlap by 5% (75%, 70%, etc.) until the endpoint of the subsequent time segment ends later than the end of the previous one.In the end, we obtained in total 1343 overlapping time segments, hence spectra, from the 51 bursts for the first stage of our spectroscopy investigation.

SPECTRAL ANALYSES & RESULTS
We performed our spectral analysis with XSPEC (version 12.12.1)using Castor Statistics (Cstat;Cash 1979).We also generated Detector Response Matrices (DRM) with the GBM Response Generator 3 released by Fermi -GBM team.As mentioned above, we fit the spectrum of each time segment (in 8−200 keV) with three different photon models: an exponentially cutoff power law model (COMPT 4 ), the sum of two black body functions (BB+BB), and a modified black body with resonant cyclotron scattering (MBB-RCS 5 ).
Based on the C-stat value obtained with each model fit, we employed the Bayesian Information Criterion (BIC; Schwarz 1978;Liddle 2007)  preferred model that is with significantly lower BIC (∆BIC > 10, which corresponds to the Bayes factor of ∼150, indicating the confidence level >99% for the likelihood ratio; Kass & Raftery 1995).
If the difference in BIC values is small, i.e., ∆BIC ≤ 10, then both models are equally preferred.
Following the fits, we also calculated the photon and energy flux of each time segment in the energy range of 8−200 keV based on the fit parameters of each photon model.Note that all errors reported throughout the paper are at the confidence level of 1σ.We also require model parameters to be well-constrained with their 1σ errors (i.e., model parameter ± 1σ error must be viable values) for the fits to be considered acceptable.

Overlapping Time-Resolved Spectroscopy
Out of the 1343 spectra modeled, we found that 1322 of them (i.e., 98.4% of the sample) are described well with COMPT, meaning that they were deemed preferred based on the BIC values.The thermal models, on the other hand, perform nearly equally in fitting: BB+BB and MBB-RCS models are statistically preferred for 542 (40.4%) and 551 (41%) out of the 1343 time segments, respectively.Here, "preferred" means that the BIC value of a model is either significantly lower than the other two models (∆BIC > 10, hence the model is the only preferred model) or comparable to the other model(s) (∆BIC ≤ 10, hence two or all three models are comparably preferred).
As stated before, our aim is to perform a clustering-based time-resolved spectral analysis to clearly reveal spectral evolution throughout short magnetar bursts.Therefore, we identified the significant spectral change points of the bursts using a machine learning based clustering algorithm with the results of spectral analysis of overlapping time segments.To do that, we chose the E peak parameter as our reference to determine the spectral change points since the COMPT model fits more than 98% of the spectra and another parameter in the model (i.e., the power-law index) does not vary significantly within individual events.
For obtaining significant change points in the sequential E peak domain, we employed the k-means clustering method (Pedregosa et al. 2011) using the midpoints of the overlapping time segments and their corresponding E peak values.With the clustering, we were able to combine the consecutive time segments that yielded similar E peak values considering their errors.Thus, we increased the statistical reliability of the spectral parameters in the second stage of spectroscopy and emphasized the E peak evolution of the burst explicitly by bringing the spectral change points to the fore.In Figure 2, we present an example of how overlapping time segments are grouped with k-means clustering for the clustering-based time-resolved spectral analysis, which reveals the spectral evolution throughout the burst, independent of initial time binning.The details of how we applied the k-means clustering to our data can be found in Appendix B.
After the k-means clustering, the grouped time segments still overlap since they are clusters of overlapping time segments (see how the last colored time segment of a cluster overlaps with the first colored time segment of the following cluster in Figure 2).Therefore, we defined the end of the cluster time interval to be the time at which half of the counts within the overlapping interval were accumulated (the black data points in Figure 2 show the non-overlapping time intervals).After finalizing the cluster time intervals, we checked the total photon counts in each time interval before subjecting them to the second-stage analysis; We found that when there are sharp changes in the spectral parameters, a non-overlapping cluster can only contain 900 (or even less in a few cases) counts.Since we confirmed that this was still statistically sufficient for the second stage of our spectral investigations, we included all spectra of non-overlapping time segments with at least 900 counts in our study.In the few cases where the burst counts remained below this level, we combined the time segments if they are adjacent, and if they are isolated single segments, we combined them with an adjacent segment with a closer E peak value.
In the end, we obtained 287 spectrally distinctive time segments for the clustering-based timeresolved spectral investigations of 51 bursts.In Table 1, we list the times of these 51 bursts, along with the detectors used in spectral analysis and event duration, as well as the number of overlapping and non-overlapping time segments.

Clustering-Based Time-Resolved Spectroscopy
We performed a detailed analysis of the 287 spectra accumulated from distinct time segments with the three continuum models (COMPT, BB+BB, and MBB-RCS).Based on these results, we determined the preferred model(s) for each segment using their BIC values.We found that COMPT is a preferred model for 279 spectra (See Table 2), out of which COMPT is the only preferred model for 83 spectra, COMPT & BB+BB are equally preferred for 98 spectra, and COMPT & MBB-RCS are equally preferred for 58 spectra.As for the thermal models, BB+BB model is statistically preferred for 145 (51%) and MBB-RCS model for 98 (34%) out of 286 spectra (one spectrum is excluded due to unacceptable fit statistics; see below about the goodness of fit test).Finally, all three models are equally preferred for 40 spectra (14%).To demonstrate the spectral shapes of these three continuum models, we present an exemplary count spectrum that nearly equally fits well with all three models in Appendix C (Figure 9).Out of the seven spectra that favor only BB+BB, we found ∆BIC ∼ 30 − 40 for four of them (two segments each from two bursts), meaning that the BB+BB model is definitely preferred over the other two models.Interestingly, the remaining time segments of these two bursts also favor thermal models (i.e., BB+BB and/or MBB-RCS), besides COMPT.In Figure 3, we present the spectral evolution of thermal model parameters from one of these two bursts.We also present the time evolution of photon flux distribution for the first seven segments of this event in Appendix C (Figure 10).To ensure that these preferred models provide statistically acceptable fits to the spectra, we also evaluated the goodness of fits for the preferred models using C-stat, following Kaastra (2017).In doing so, we used the C-stat values excluding 30−40 keV to avoid the contribution from the iodine K-edge, which could affect the statistics for bright events6 .To this end, we computed7 the expected C-stat (C e ) and its variance (C v ) based on model predictions and determined whether each model fit remained within the 3σ level8 of its expected C-stat distribution.We found that for only one spectrum (the second time segment of the burst at 608204639.277MET), all three models yielded C-stat values with > 3σ deviation; this is due to the counts in < 10 keV being lower than what is expected from these three models.Therefore, we excluded this time segment from our analysis and evaluated 286 out of 287 time segments from the 51 bursts in what follows below.Note that the numbers given in Table 2 exclude this time segment.
In Figure 4 panel (a), we present the scatter plot of the 279 pairs of photon index (Γ) and E peak values of the COMPT model.Note that the plot is color-coded based on the energy flux in 8−200 keV.We find that the distribution of Γ is asymmetric, described best by a Gaussian with an underlying first-order polynomial (see panel (b) of Figure 4).Γ values above 0.25 follow a Gaussian with a mean value of 0.65±0.02and σ = 0.18±0.02.However, Γ values below 0.25 form an excess above the Gaussian tail, which can be described with a first-order polynomial of (7.08±1.95)Γ+ (10.13±1.76).It is clear from Figure 4 that those Γ values forming the excess are obtained from time segments with the lowest flux in our sample (below 1 × 10 −5 erg cm −2 s −1 ).These time segments with low flux values generally correspond to the first or last time segments of the bursts as expected.On the other hand, time segments with higher flux values yield Γ values mostly larger than 0.25, and those are the values whose distribution is consistent with a Gaussian.It is also important to note that the time segments with the highest flux values (above 4 × 10 −5 erg cm −2 s −1 ) yield Γ values in a very narrow range from 0.2 to 0.55.As for the corresponding E peak values, they range from ∼20 to ∼52 keV, and their distribution follows a Gaussian shape with a mean of 34.39±0.26keV and a width, σ = 4.21±0.20 (see panel (c) of Figure 4).Moreover, there appears a positive correlation between E peak and flux.To quantify the correlation, we computed Spearman's rank order correlation coefficients (ρ and chance probability, P ) for E peak and flux using a bootstrap method; we generated 10000 data sets by taking into account the uncertainties of E peak and flux, and calculated the correlation coefficients for each data set.From the distribution of these coefficients, we obtained the mean and 1σ confidence interval of ρ and P .
We found ρ = 0.81 ± 0.01 and P < 10 −63 for E peak and flux, which lends support to the observed positive correlation.
For the 145 BB+BB-preferred spectra, we present the scatter plot of low BB (kT Low ) and high BB (kT High ) temperatures in Figure 5.It is again color-coded based on the energy flux in the 8−200 keV.We find a positive correlation between the two parameters (ρ = 0.67 ± 0.04 and P < 10 −17 ).Overall, lower temperatures are associated with low flux levels (purple and blue data points in panel (a) of Figure 5), while higher temperatures correspond to high flux values (red and yellow data points).The intermediate flux values, however, have a much wider BB+BB temperature range.In terms of parameter distributions, we find that a best Gaussian fit to the distribution of kT Low yields a mean of 6.30±0.17keV with a sigma of 1.80±0.15(see panel (c) of Figure 5).On the other hand, our criterion that kT High parameter must be larger than kT Low parameter forces the distribution of kT High to be asymmetrical (see panel (b) of Figure 5).Therefore, the best fit to the distribution of kT High results in a two-sided Gaussian fit.The distribution has a mean of 8.81±2.74keV with σ left = 0.2±1.8 and σ right = 4.83±1.4.Finally, for the 98 MBB-RCS-preferred spectra, the distribution of temperature (kT M ) lies between 5−13 keV, which remains in between kT Low and kT High distributions of the BB+BB model.In Figure 6, we present the kT M parameter distribution of the MBB-RCS model.The distribution is described well with the normal distribution, which peaks at 7.84±0.12keV with σ = 1.12±0.1.As in the BB+BB model, we also observe a positive correlation between energy flux and temperature in this model (ρ = 0.62 ± 0.02 and P < 10 −11 ); on average, higher kT M values correspond to higher flux values.

DISCUSSION
In this study, we applied a novel approach for the first time in magnetar bursts research: a clusteringbased time-resolved spectroscopy.We performed high time-resolved spectral analysis of the brightest SGR J1935+2154 bursts by first dividing them into the sequentially-overlapping shortest time intervals to clearly reveal the spectral evolution during the short magnetar bursts.We then identified the significant spectral change points throughout the bursts using the spectral analysis results of these overlapping time segments through a machine learning based clustering approach and completed the second round of spectroscopy with the data extracted from the time intervals between these change points.In the end, we obtained 287 non-overlapping sequential time segments from 51 bursts.Out of 287 spectrally-distinguishing time intervals, 207 arise from the brightest half of our burst sample, indicating the existence of a quite significant spectral evolution in the brightest magnetar bursts.Lin et al. (2020a,b) studied time-integrated spectral properties of 275 SGR J1935+2154 bursts between 2014 and 2020 observed with Fermi -GBM.They found that ∼62% and ∼38% could be described well with the BB+BB and COMPT models, respectively.On the contrary, our clustering-based time-resolved spectral analysis revealed that only ∼51% can be fit with BB+BB while nearly all time segments (∼98%) of the brightest bursts could be represented with COMPT.Note the fact that our burst sample selected from these 275 bursts are the brightest ones, based on our criterion as explained in Section 2; therefore, it would be a fairer comparison with our results if we exclude dim bursts from their statistics.In that case, the accepted fit percentages for their samples (178 bursts) increase to ∼96% and ∼58% with the BB+BB and COMPT models, respectively, which still are quite different from our results found here.Besides the difference between time-integrated and time-resolved analysis, such a difference in model preferences of the two studies might arise from the fact that our sample includes also SGR J1935+2154 bursts from 2021-2022 active episode.However, our statistics remained nearly the same when we excluded these bursts and evaluated only the bursts before 2021, and even when we look at the percentages of statistically acceptable fits (instead of BICbased preferences), COMPT is a better-describing model than BB+BB.We thus conclude that the time-resolved spectral results are different from those of time-integrated spectral investigations.This is not surprising since we observe significant spectral variations throughout bursts and superposition of spectra with varying E peak of COMPT could mimic the spectral energy distribution of BB+BB or results in the distribution that deviates from a single-component COMPT model.Besides, the timeintegrated spectra naturally have more counts even at higher energies, with which the two-component model parameters are better determined.

Spectral Parameters-Flux-Area Correlations
As for the model parameters of the time-integrated spectroscopy of SGR J1935+2154 bursts between 2014−2016, those that fit well with the COMPT model have E peak values range between ∼25 and 40 keV with a Gaussian mean of 30.4±0.2 keV (Lin et al. 2020a).The later 2019−2020 bursts have a wider range of E peak (∼10−40 keV) with a slightly lower mean of 26.4±0.6 keV (Lin et al. 2020b).In comparison, our clustering-based time-resolved spectroscopy yields a slightly higher energy range of ∼20−52 keV and a higher mean E peak of 34.4±0.3keV.We again note here that our burst sample includes bursts detected in the 2021 and 2022 active episodes and met our brightness criterion; when we compare only the bright events between 2016 and 2020, the time-resolved E peak in our sample is still harder.
Previously, time-integrated spectral analysis of SGR J1935+2154 bursts and other prolific magnetar bursts (SGR J1550−5418) showed that E peak is correlated with the energy flux or fluence described by a power law or a broken power law (van der Horst et al. 2012;Lin et al. 2020a,b).We present in Figure 7 (left panel), the E peak vs. flux plot for our time-resolved results; the correlation was revealed more clearly as a result of our time-resolved spectroscopy (Spearman's rank correlation coefficient, ρ = 0.81 ± 0.01, P < 10 −63 ).To better quantify the relation between E peak and flux, we fit the trend as follows: Since flux errors are small (∆F/F ≲ 0.015), we grouped the data in the flux domain such that each group would include 20 data points.For each group, we computed the weighted mean flux and E peak , as well as 1σ uncertainty of E peak .Modeling the grouped trend with a single power law model (PL) yields an unacceptable fit (χ 2 /dof = 117.4/12).A fit with a broken power law model (BPL), on the other hand, results in statistically acceptable representation (χ 2 /dof = 15.7/10),yielding a break at the flux of (2.00±0.05)×10 −5 erg cm −2 s −1 , and positive indices of 0.08±0.01 and 0.37±0.03before and after the break, respectively (see Figure 7 left panel).
On the other hand, we found that the range of the COMPT photon index (Γ) parameter of our timeresolved investigation is consistent with the time-integrated spectroscopy although their distributions are quite different.The Γ distribution of time-integrated spectral analysis runs from −1.5 to 1 and follows a Gaussian with a mean of ∼0 (Lin et al. 2020a,b).In the case of our time-resolved spectroscopy, the distribution has a tail, best described with a Gaussian with an underlying firstorder polynomial.The positive Γ values above 0.25 are distributed like a Gaussian with a mean at 0.65±0.02,while Γ values below 0.25 form an excess above the Gaussian tail.We found that these indices (Γ < 0.25) were obtained from the spectra with the lowest flux values in our sample (< 1×10 −5 erg cm −2 s −1 ).Moreover, the spectra with the highest flux values in our sample (> 4 × 10 −5 erg cm −2 s −1 ) yield Γ values in a very narrow range with a weighted mean of 0.47±0.03.There exists a positive correlation between Γ and flux up to a certain flux level that coincides with the flux break of E peak vs. flux (see the right panel of Figure 7).The index and flux are anti-correlated after that point.The change in trend between index and flux is also indicated in Figure 4: the photon index increases with increasing flux up to about 2 × 10 −5 erg cm −2 s −1 (blue and purple data points in the panel (a) of Figure 4).Then, the index starts to decrease with increasing flux.In modeling this trend, we followed a similar approach as in the case of E peak vs. flux: We obtained weighted mean Γ values and their 1σ uncertainties for the groups of 20 flux values and used a broken log-linear function to fit, yielding χ 2 /dof = 10.1/10.We find the break at the flux of (2.04±0.09)×10−5 erg cm −2 s −1 , that is consistent with the E peak vs. flux case.The slope of 0.97±0.09before the break changes to −0.86±0.17afterwards (see the right panel of Figure 7).Younes et al. (2014) performed time-resolved spectral investigations of bursts from another prolific magnetar, SGR J1550−5418, also observed with Fermi -GBM.They also found that the E peak vs. flux relation is described better with a BPL rather than a single PL.This break point is at the flux of ∼ 1 × 10 −5 erg cm −2 s −1 , which is about a half the break value in flux that we found for SGR J1935+2154.Unlike our findings, the E peak vs. flux correlation of SGR J1550−5418 is negative at low flux values, while it is positive after the flux break.Moreover, Γ remains constant at ∼ −0.8 up to the flux break, then follows the same positive trend as seen between E peak and flux.Similar dual relation between the E peak and flux with a break at the flux of ∼ 1 × 10 −5 erg cm −2 s −1 was also reported in the time-resolved spectroscopy of five bright bursts from yet another magnetar, SGR J0501+4516 (Lin et al. 2011).
For the BB+BB model, earlier time-integrated spectral studies of SGR J1935+2154 bursts yielded parameters of ∼2−8 keV for kT Low with a mean of 4.5 keV, and ∼8−20 keV for kT High with a mean of 11 keV (Lin et al. 2020a,b).Our clustering-based time-resolved spectroscopy using the BB+BB model results in similar parameter range and distribution for kT High .However, our kT Low parameter reaches up to 10 keV with a larger mean of 6.3±0.2 keV.Moreover, unlike the time-integrated spectral analysis, our study reveals a positive correlation between kT Low and kT High (ρ = 0.67 ± 0.04, P < 10 −17 ; see panel (a) of Figure 5).
Based on the BB+BB parameters, we also calculated the size of blackbody emitting areas as R 2 = (F d 2 )/(σT 4 ) in km 2 where F is the flux, d is the distance to the source (here we use a distance of 9 kpc, consistent with Lin et al. 2020a,b), σ is the Stefan-Boltzmann constant, and T is the temperature in Kelvin.Our study revealed that the relation between T and R 2 varies with the burst flux, similar to the findings of time-resolved spectral analysis of SGR J1550−5418 bursts (Younes et al. 2014).Therefore, following Younes et al. (2014), we divided our 145 spectra that favors BB+BB model into four groups based on energy flux (in erg cm −2 s −1 ): F < 10 −5.5 , 10 −5.5 < F < 10 −5.0 , 10 −5.0 < F < 10 −4.5 , and F > 10 −4.5 .In the left panel of Figure 8, we present the plot of R 2 vs. kT for the above-mentioned flux ranges.We found that the best fit for the R 2 vs. kT relation of the three highest-flux data groups is the BPL model.However, for the lowest-flux group, both PL and BPL provide statistically acceptable fits.Despite the difference, we observe that the overall trends in all four groups are the same: The negative correlation (slope) between R 2 and kT becomes steeper after the break of the BPLs.We summarize our flux-dependent R 2 vs. kT fit results in Table 3.Note that these fit results are obtained with the data from all of the spectra that favor BB+BB, not with the weighted means of the data (that we show in Figure 8 only for display purposes).
We observe that the kT range of SGR J1935+2154 and SGR J1550−5418 bursts are similar while R 2 for SGR J1935+2154 bursts do not go down below 1 km 2 unlike SGR J1550−5418 bursts.Still, the relation between R 2 and kT has the same trend for both sources: R 2 decreasing more slowly for kT Low compared to kT High across all flux ranges.We found that the slope between R 2 & kT Low (α-kT Low ) and the slope between R 2 & kT High (α-kT High ) for all flux groups in both studies are consistent with one another within at most 3σ errors.Also, kT break values are consistent with one another within errors (∼7−9 keV for SGR J1550−5418; Younes et al. 2014).Moreover, we observed a similar relation between R 2 and kT Low : α-kT Low becomes steeper as flux decreases.However, our investigations do not yield a decreasing trend in the α-kT High with decreasing flux.Furthermore, similar to SGR J1550−5418 bursts, a single PL model yields statistically acceptable R 2 & kT fit in the lowest flux group.
Finally, for the other thermal model, namely, the modified blackbody whose emission undergoes resonant cyclotron scattering (MBB-RCS), we observe a correlation between temperature and flux Table 3. R 2 vs. kT fit parameters (PL index α and break energy) for various flux ranges of BB+BB and MBB-RCS models as shown in Figure 8 BB+BB MBB-RCS −1.02 ± 1.64 −4.88 ± 0.79 10.08 ± 1.12 −2.01 ± 0.10 a 10 −5.0 < F < 10 −4.5 −1.14 ± 0.35 −5.62 ± 0.24 8.98 ± 1.02 −2.01 ± 0.10 a 10 −5.5 < F < 10 −5.0 −1.36 ± 0.81 −4.99 ± 0.47 6.96 ± 1.09 −3.83 ± 0.15 F < 10 (see Figure 6, ρ = 0.62 ± 0.02 and P < 10 −11 ).Also for the 98 MBB-RCS-preferred spectra, we calculated the corresponding emitting region (R 2 ) based on their temperature (kT M ).In the right panel of Figure 8, we present the scatter plot of R 2 vs. kT M with color-coding based on the same flux ranges used for BB+BB − except that we combined the highest two flux groups for the MBB-RCS model due to the deficiency of time segments (only three) in the highest flux regime well described by MBB-RCS.We find that for R 2 vs. kT M the best fit across all flux ranges is a PL, which is presented also in Table 3.We also find that PL indices (α-kT M ) for the lowest two flux regimes are consistent with −4 (i.e., as expected from F ∝ σT 4 ).More importantly, α-kT M obtained from the highest flux spectra (F > 10 −5 erg cm −2 s −1 , corresponding to an isotropic luminosity of 10 41 erg s −1 ) is significantly different, −2.01 ± 0.10, than the lower-flux spectra.

Interpretative Elements
We now elaborate on the interpretation of some of our results.Large Thomson opacities are expected for magnetars in outburst, quickly discernible using non-magnetic estimates that were identified by Baring & Harding (2007).Setting E e ≳ L γ /(4πR 2 c) as the representative kinetic energy density in radiating electrons, then if they possess a typical Lorentz factor ⟨γ e ⟩ ∼ 1 , one quickly arrives at an electron number density n e ∼ E e /(m e c 2 ) , so that the non-magnetic Thomson optical depth is τ T = n e σ T R ≳ L γ σ T /(4πRm e c 3 ) ; this is the familiar compactness parameter.For R ∼ 10 7 cm, this yields τ T ∼ 10 4 (i.e., n e ≳ 10 21 cm −3 ) for SGR bursts of typical isotropic luminosities L γ ∼ 10 40 erg s −1 , indicating optically thick, super-Eddington conditions (e.g.Thompson & Duncan 1996) that drive plasma flow along the field lines.
In Figure 8 (left), we observe a significant deviation from the Stefan-Boltzmann (S-B) law for an isotropic radiation field (R 2 ∝ kT −4 ): We find R 2 ∝ kT −α , where α ∼1−1.4 for kT ≲ 7 − 10 keV above the flux level of ∼ 10 −5.5 erg cm −2 s −1 .Note that this flux corresponds to an isotropic luminosity of 3 × 10 40 erg s −1 at a source distance of 9 kpc.The resulting broken power-law R 2 − kT correlations can be interpreted as a signature of the spatial extension of the active emission region, which can be presumed to be a broad, flaring flux tube.If this tube has a transverse dimension of R t ∼ 2 − 10 km for its cross section at the highest altitudes, one can infer a tube length of R l ∼ 4R 2 /R t ∼ 100 − 500 km for the largest R 2 values.The highly optically thick gas in the emission region will naturally cool adiabatically when moving between smaller, hotter regions near the flux tube footpoints at the stellar surface, and the high altitude, large area regions near the equatorial tube apex.This yields the spectral extension we see.If the radiating gas is locally quasi-thermal, then near the footpoints the magnetic field is high, and photospheric/atmospheric simulations (see Fig. 2 of Hu et al. 2022) of radiative transfer in this sub-cyclotronic frequency domain (ω ≪ ω B = eB/ℏc) indicate that the radiation is quasi-isotropic.In contrast, since the magnetic field is much lower (likely 3-4 orders of magnitude) near the tube apex, the Compton scattering radiative transfer in the local region samples the cyclotronic domain where ω ∼ ω B .This domain evinces significantly anisotropic emergent radiation fields (see the right column of Fig. 2 of Hu et al. 2022), with a decrement of intensity over a wide range of directions that are oriented closer to the outer surface of the magnetic flux tube/photosphere; this tube surface is aligned with the local field direction that guides plasma motion and associated adiabatic cooling.Such a decrement alters the S-B law from its isotropic 'R 2 (kT ) 4 = constant' form, and lowers the perceived area in lower kT regions at higher altitudes9 .The reduced values of R 2 (kT ) 4 naturally generate the breaks apparent in Figure 8.In particular, the actual flux tube areas would rise above those the S-B estimate obtains.Concomitantly, this decrement yields the break in the flux-E peak correlation depicted in Figure 7 (left).
The extended spatial emission zone scenario suggests that a multi-blackbody fit would be preferable to just a BB+BB one.This may be the reason why the MBB-RCS model provides a good spectral fit for a significant portion of the sample studied here, in which for the first time the MBB-RCS model has been applied to a substantial sample size.Clearly, the scattering component of the MBB-RCS picture applies outside the highly optically thick primary burst emission zone.If the plasma density becomes high enough, as may be more likely for the highly energetic events (L ≳ 10 41 erg s −1 ), multiple scatterings of photons by the magnetospheric charges would arise, and this likely would harden the emerging spectrum.The required plasma densities would be considerably higher than those needed in resonant inverse Compton emission models (Baring & Harding 2007;Fernández & Thompson 2007;Wadiasingh et al. 2018) of the persistent hard X-ray tail emission of luminosities L ∼ 10 35 erg s −1 (Kuiper et al. 2004;Götz et al. 2006;den Hartog et al. 2008) from various magnetars.In such domains, the MBB-RCS model would need to be expanded (Yamasaki et al. 2020).
Accordingly, strong motivations exist for future detailed accounting of the altitudinal dependence of the anisotropy of the flux tube radiation field, in combination with RCS modeling in neighboring magnetospheric regions.The RCS process would have to address both single and multiple scattering domains.This level of sophistication is desirable for more precisely interpreting area-color correlations at different luminosity levels.In particular, if the energies kT ∼ 7 − 10 keV of the breaks are interpreted as a loose measure of the ℏω B value (when the anisotropy becomes substantial) somewhat near the tube apex, it may prove possible to constrain the active flux tube dimensions and magnetospheric locale using the area-color correlations.Thus, our results not only encourage further employment of the MBB-RCS physical model in diverse burst samples from various sources but also highlight the need for more comprehensive modeling approaches in understanding the behavior of highly energetic magnetar flares.
A. SPECTRAL DATA EXTRACTION For each event, we calculated the burst duration (i.e., the duration over which burst spectral extraction is to be performed) using the data collected with the brightest NaI detector, which is the detector with the smallest zenith angle to the sky position of the source.In particular, we first constructed a light curve with 4 ms time resolution for the time interval of −10 s to 10 s with respect to the burst start time published in Lin et al. (2020a,b), using time-tagged photon data in the 8−200 keV energy band.We then generated the Bayesian Block representation of the light curve.Blocks longer than 4 s were considered background blocks, and the background level was calculated by averaging the count rates of these background blocks.The blocks shorter than 4 s and higher in rates than the background level are taken as burst blocks, and the duration is calculated as the time interval from the start of the first burst block to the end of the last one (Lin et al. 2013).
For defining sequentially overlapping time segments of each burst, we first determined the background level as the average of the pre-burst time interval between −50 s and −1 s by taking the start time of Bayesian Block duration as the reference point.Then, we obtained a background-subtracted light curve with 4 ms time resolution in the energy range of 8−200 keV and determined the beginning and end points of time segments for the spectral data extraction, each of which contains at least 1200 background-subtracted counts and overlaps 80% in time with the previous one.As stated in Section 2, in the case of accumulation of time segments at the peak of a burst, we reduced the overlap by 5% incrementally until the endpoint of the subsequent time segment ends later than the end of the previous one.
We aim to investigate each burst by dividing it into as many short time intervals as possible.At the same time, we also aim for these intervals to have sufficient statistics (i.e., burst counts) for reliable spectral analysis.Therefore, setting a threshold background-subtracted counts to be included in each time segment is required.To this end, we selected a small set of bursts and extracted spectra of the time segments for each of these bursts, again with 80% overlap, but with a set of different numbers of background-subtracted counts, namely 600, 800, 1000, 1200, 1500, 1800, and 2000.We then fit these spectra with the three models that we used in this study and checked whether the resulting model parameters are constrained within at least 2σ level.Based on this analysis, we concluded that a threshold of 1200 background-subtracted counts was needed to ensure reliable spectral analysis results with constrained model parameters for each time segment.With this requirement of 1200 burst counts in each time segment, the minimum number of counts throughout a burst should be larger than about 2400 in order for a burst to consist of at least two non-overlapping time segments, therefore, spectral evolution could be defined.
We also checked the potential count saturation in the data, from which very bright bursts may suffer.The GBM TTE data are affected by 2.6 µs deadtime as a result of the fixed data-packet processing speed of 375 kHz on the spacecraft (Meegan et al. 2009).Therefore, when the combined count rates of all GBM detectors surpass this limit, data saturation occurs.We did not find saturation in the data during any of the 51 bursts in our chosen sample.

B. MACHINE LEARNING APPROACH FOR SPECTRAL CLUSTERING
K-means clustering is a commonly used algorithm in data analysis that groups data points based on specific features, which in our case are the midpoints of time segments and corresponding E peak values for each burst.Its primary objective is to cluster data points in a way that maximizes the similarity within clusters and minimizes the similarity between different clusters (Lloyd 1982).
In this study, Python programming language (version 3.6.9)and Scikit-learn (version 1.2.1) were used for k-means clustering.The algorithm requires a predetermined number of clusters (k).Firstly, it initializes centroids by randomly selecting k samples from the data set.After initialization, the algorithm iterates between the remaining two steps: It first assigns each sample to the nearest centroid, then, creates new centroids by computing the mean value of all the samples assigned to each of the previous centroids.The algorithm repeats these last two steps until the squared difference between the old and new centroids is less than a predefined threshold, indicating that the centroids have become stable, and the clustering is complete.
For each burst, first, the data was scaled to prevent the effect of one variable from overriding the other since the ranges of time and E peak axes vary significantly.Next, the k-means algorithm was implemented for all possible k values ranging from one to (N − 1), where N is the number of time segments, and the corresponding inertia (i.e., the sum of squared distances of samples to their closest cluster center) for each k was recorded.Note that the reciprocal of E peak errors was given as a weight for the inertia calculation.By definition of inertia, its value decreases as the number of clusters (k) increases.However, increasing k does not significantly reduce inertia after a certain number of k, and the inertia vs. k graph becomes flat after this point, encompassing a minimum of about 20−25% of the high-k values in our sample.After studying a sample of events, we found that this point corresponds to the optimal number of clusters for finding the largest spectral variations during each burst.Therefore, instead of heuristic techniques that are commonly used to find the optimal k in k-means clustering, e.g., the "elbow" method, we developed our method as follows: First, we calculated the average of inertia values that corresponds to the highest 25% of the k values in the inertia vs. k graph.We then added 1% of the maximum inertia value to the average inertia of the highest 25% of the k values.The nearest integer k value corresponding to the resulting inertia is then used for the optimal number of clusters.
Finally, we note that in addition to k-means clustering we also tested other machine learning based clustering algorithms; namely, density-based spatial clustering, agglomerative clustering, and Gaussian mixture model.They all require a number of clusters to be specified and their results are consistent with those of k-means.We, therefore, consider k-means a robust method for our purpose.

C. CONTINUUM MODEL COMPARISONS
To demonstrate how different photon models represent observed burst spectra, we present in Figure 9 the spectrum extracted from a time segment of the burst that occurred at 488642074.718(Fermi MET; 160626 13:54:30.722UTC).All three models employed in this investigation yield equally well fit to this particular burst data in our energy passband: MBB-RCS with kT of 9.46 keV (red curve in Figure 9), COMPT with Γ = 0.56 and E peak = 39.07 keV (blue curve), and BB+BB with kT Low = 7.4 & kT High = 15.31keV (green curve).Slight differences among the models are observed at the low-energy end and at the high-energy end, which are statistically indistinguishable.
We also present the evolution of spectral parameters of another example burst detected at 652927551.870(Fermi MET; the same event also shown in Figure 2 and 3) in Figure 10.This is an example event, throughout which a thermal model is preferred.In the left panel, we present the kT parameters of the thermal models (i.e., BB+BB and MBB-RCS) for the first 7 time segments.In the right panel, we present the evolution of photon flux distribution with time (from yellow to red) using the BB+BB model for segments 2 through 7. Note that both COMPT and MBB-RCS model curves are displayed for the first time segment.Since photon flux distributions of segments 2 to 7 are similar at high energies (above 60 keV), we zoomed into the energy range of 8 − 60 keV to clearly exhibit model differences at low energies.
Figure1.Light curve of an SGR J1935+2154 burst detected at 2021 January 29 10:35:39.918UTC as seen with the brightest detector (n5).The vertical dashed lines show Bayesian Block duration start and end times, respectively.The red horizontal lines represent the 22 overlapping time segments with each subsequent segment having an overlap of 80%.We note that the vertical values corresponding to the red lines are arbitrary, and chosen only for display purposes.

Figure 2 .
Figure 2. Light curve of an SGR J1935+2154 burst detected on 2021 September 10 at 00:45:46.875UTC as seen with the brightest detector (n8) is shown with grey dotted lines (right axis).Via k-means clustering, 160 overlapping time segments and their corresponding E peak values (the colored data points with asymmetric E peak error bars) yielded 13 spectrally distinctive clusters, each of which is shown with a different color.The black data points show the E peak values with asymmetric error bars that are obtained from the COMPT fit of the data extracted from the time spanned by these 13 clusters in the second stage of spectral analysis.

Figure 3 .
Figure3.Blackbody temperature evolution for a burst, throughout which a thermal model is preferred.The light curve is shown with grey dotted lines (right axis; the same event as Figure2).Thick data points show the kT Low and kT High parameters of BB+BB while thin data points with circles show kT M of MBB-RCS.The color code represents the preferred model(s) for each time segment based on ∆BIC.The ∆BIC between BB+BB and the other two models is only slightly above 10 for the first time segment, which is shown with blue color.

Figure 5 .
Figure 5. (a) The scatter plot of kT Low vs. kT High parameters of 145 spectra that can be described well with BB+BB.The colors indicate the energy flux values.(b) The distribution of kT High with the best fit model (two-sided Gaussian) curve is shown in red.(c) The distribution of kT Low values with the best fit Gaussian function overlaid in red.

Figure 6 .
Figure 6.The distribution of kT M parameter of 98 spectra that favor MBB-RCS, with the best-fit Gaussian function overlaid in red (left axis).The scatter plot of energy flux vs. kT M is also shown (right axis).The colors indicate the energy flux values.

Figure 8
Figure 8. [Left] Flux color-coded plot of R 2 vs. kT for 145 BB+BB spectra.Each data point represents the weighted means of R 2 and kT of three time segments only for display purposes.Solid lines show BPL fits.The dark blue dash-dotted line indicates PL fit to the lowest flux group.The lowest flux group is fitted almost equally well with both a single PL and BPL.The black dashed line indicates R 2 ∝ kT −4 .[Right] Flux color coded scatter plot of R 2 vs. kT for 98 spectra favoring MBB-RCS.Solid lines represent PL fits.The dashed line indicates R 2 ∝ kT −4 .Note that we combined the data points of the two highest flux groups for MBB-RCS since only three spectra in the highest flux group (F > 10 −4.5 erg cm −2 s −1 ) favor the MBB-RCS model, which are shown with red diamonds.
Slopes of the data points from the highest two flux groups for the MBB-RCS model are the same since they are combined due to the deficiency of time segments (only three) that favor MBB-RCS in the highest flux regime.b A single PL fit to the data.

Figure 9 .
Figure 9. Observed count spectrum of a time segment from the burst detected at 488642074.718(Fermi MET) are shown in black crosses.The three solid trends are the best-fitting model curves: MBB-RCS in red, COMPT in blue, and BB+BB in green, all of which fit the spectrum equally well.
Yamasaki et al. (2020)odel to be used in XSPEC by following the prescription ofYamasaki et al. (2020).The table model covers the energy range from 5 to 300 keV, and the parameter grid of T eff consists of 79 values between 1 and 40 keV with increments of 0.5 keV.This table model can be downloaded from https://zenodo.org/records/10485159.

Table 1 .
List of 51 SGR J1935+2154 bursts included in our sample The bursts in 2016, in 2019−2020, and in 2021−2022 are taken from Lin et al. (2020a), Lin et al. (2020b), and Lin et al. (in preparation), respectively.b Unblocked NaI detectors used in spectral analysis.The brightest detectors shown in bold are used to determine the start and stop times of time segments for the extraction of spectra.c Bayesian Block Duration.Event times (both UTC and MET) represent the Bayesian Block duration start times of the bursts.
Table 1 continued on next page Table 1 (continued)

Table 2 .
Number of non-overlapping time segments per preferred photon models Note-a Number of model parameters is indicated in parentheses.