Fermi/GAMMA-RAY BURST MONITOR OBSERVATIONS OF SGR J0501+4516 BURSTS

For the T₉₀¹⁸ estimate we used the algorithm originally developed by Kouveliotou et al. (1993) for gamma-ray bursts (GRBs) and later adapted for SGR bursts (Göğüş et al. 2001; Gavriil et al. 2004), modified slightly to accommodate the GBM TTE data. Each duration was calculated between 8 and 100 keV with 2 ms bins as follows. First, we fit the burst background using two time intervals before and after each burst ([−2 s, −0.5 s] and [0.5 s, 2 s]) with a first-order polynomial; these intervals were kept mostly the same for all bursts unless a precursor and/or a tail were detected. Then we fit the background-subtracted, cumulative burst counts, with a linear plus a step function simultaneously; the linear part was user-selected (before and after the burst), while the step part was determined from the fit linear trend. The height of the step function was then used to represent the net total counts ($\textit {N}$) of the burst, subsequently used for the T₉₀ determination.

Panels (a) and (b) in Figure 2 show the distributions of T₉₀ and T₅₀ for all 29 bursts; all individual T₉₀ values can be found in Table 1 (Column 5). For comparison with other magnetar duration distributions (see Section 3.4) we fit each distribution with log-normal functions and obtained 〈T₉₀〉 = 122.6^+7.9_{− 7.5} ms (σ = 0.35 ± 0.03, where σ is the width of the distribution in the log frame) and 〈T₅₀〉 = 31.6^+2.5_{− 2.3} ms (σ = 0.30 ± 0.03). The average values of the raw data weighted by their errors are 〈T^w₉₀〉 = 138.3_−20.5^+1.07 ms and 〈T^w₅₀〉 = 32.4_−0.8^+0.9 ms.

Zoom In Zoom Out Reset image size

The distributions above did not account for the uncertainties in the durations. We estimate the (max, min) errors ($\Delta {T_{\rm 90{\rm max}}}$, $\Delta {T_{\rm 90{\rm min}}}$) for each T₉₀ (similarly T₅₀) as

$$\begin{eqnarray} T_{90{\rm max}} &=& t_{0.95(N+\sqrt{N})} - t_{0.05(N-\sqrt{N})} \nonumber \\ T_{90{\rm min}} &=& t_{0.95(N-\sqrt{N})} - t_{0.05(N+\sqrt{N})} \nonumber \\ \Delta {T_{90{\rm max}}}&=&{T_{90{\rm max}}}-{T_{90}} \nonumber \\ \Delta {T_{90{\rm min}}}&=&{T_{90}}-{T_{90\rm {{\rm min}}}}, \end{eqnarray} \tag{ 1 }$$

where t_n is the time when the cumulative light curve reaches n counts. We then created probability distribution functions (pdfs) for each T₉₀ and T₅₀, as described by Starling et al. (2008) and Evans et al. (2009). Each pdf was constructed using a two-sided normal distribution, where the width of each half (σ) was set to the uncertainty of the duration ($\Delta {T_{\rm 90{\rm max}}}$, $\Delta {T_{\rm 90{\rm min}}}$, $\Delta {T_{\rm 50{\rm max}}}$, and $\Delta {T_{\rm 50{\rm min}}}$):

$$\begin{eqnarray} P(x|\bar{x}, \sigma _1, \sigma _2) = \frac{\sqrt{2}}{\sqrt{\pi }(\sigma _1 + \sigma _2)} \left\lbrace \begin{array}{@{}ll}e^A & (x \le \bar{x}) \\ e^B & (x > \bar{x})\\ \end{array}\right.\cr {\rm with}\ \left\lbrace \begin{array}{@{}ll}A={-(x - \bar{x})^2 / 2\sigma _{1}^{2}} \\\ms B={-(x - \bar{x})^2 / 2\sigma _{2}^{2}}. \end{array}\right. \end{eqnarray} \tag{ 2 }$$

Each pdf describes the likelihood to obtain x (i.e., T₉₀ or T₅₀) given its measured value $\bar{x}$. Finally, we averaged all sample pdfs to create the total pdf of T₉₀ (T₅₀) shown in panels (a) and (b) of Figure 2. Note that the shorter (fainter) events have larger errors, resulting in "pulling" the pdf toward shorter durations.

The photon-based durations, T^ph₉₀, are estimated with an algorithm similar to the one used above over each burst cumulative fluence in erg cm⁻². We used the same time resolution (2 ms) and energy range (8–100 keV) as in the count durations. The essential difference here is that these measurements utilize the intrinsic (deconvolved) burst spectra instead of the detector recorded counts to define the burst intrinsic durations independent of different instruments. To perform these estimates, we used the GBM public software tool RMFIT version 3.3¹⁹ (for a description of this tool see also Kaneko et al. 2006) and the new data type CTTE specially created to facilitate analyses of short events. This data type simply bins the 128 TTE energy channels into the same eight bins as the CTIME data. The errors in the duration estimates are taken from Koshut (1996) and Koshut et al. (1996).

A detailed description of the photon-based durations can be found in the First Two Years GRB Catalog of Fermi/GBM (W. Paciesas et al. 2011, in preparation). In short, an adequate background interval is selected before and after each burst and fit with the lowest acceptable order of a polynomial to determine the background model parameters. Next, the entire burst interval is fit to determine the default set of photon model parameters. The model used in these fits is a power law with an exponential cutoff (COMPT; described in detail in Section 4). When all background and source model selections are determined for each 2 ms time bin, we subtract the background, fit its spectrum using the COMPT model, and calculate its photon flux. These values are then used as inputs for the T^ph₉₀ (T^ph₅₀) estimates, performed with the same algorithm described above.

Figure 3 shows the distributions of Tph90 (Tph50) fit with a log-normal function (panels (a) and (b)), obtaining 〈Tph90〉 = 124.2−15.2+17.3 ms (σ = 0.38 ± 0.06, where σ is the width of the distribution in the log frame) and 〈Tph50〉 = 27.6−1.7+1.8 ms (σ = 0.21 ± 0.03). The average values of the raw data weighted by their errors are 〈Tphw90〉 = 161.2+1.6− 1.6 ms and 〈Tphw50〉 = 49.2+0.8− 0.8 ms. The individual Tph90 values can be found in Table 1 (Column 6).

Zoom In Zoom Out Reset image size

The solid curves in Figure 3 exhibit the pdfs for T^ph₉₀ and T^ph₅₀ after taking into account the symmetrized errors in the values. This plot is an adaptive kernel density estimation: it was made by adding up a set of normalized Gaussian functions, one at each of the data points and with a 1σ width given by the corresponding error estimate. The multiple narrow spikes in this figure, mostly at the long-duration end of the plot, are largely due to the fact that the corresponding errors (in contrast to those of shorter and fainter durations) are relatively small, yielding tall and narrow Gaussian components. This is a rather high-variance estimate of the distribution function. As an alternative, we show in Figure 4 a completely different density estimate based on an adaptation of the Bayesian block algorithm for histogramming of the T^ph₉₀ and T^ph₅₀ data, taking into account measurement errors in the independent variable (Scargle 1998; J. Scargle et al. 2011, in preparation). This algorithm finds the optimal piecewise constant model to represent the data; the optimization corresponds to the maximum likelihood for a constant-rate Poisson model for the data in each of the bins. The sizes and locations of the bins are all determined by this optimization, not pre-defined as in ordinary histograms. The solid curve is the Bayesian block representation superimposed on such an ordinary histogram with bins chosen so that the value plotted is more than one only if there are duplicate values. This analysis suggests that the spiky structure of the pdfs in Figure 3 is mostly due to noise fluctuations (the spikes corresponding to a small number of points with small formal measurement errors) and does not support a quantized or a multimodal duration distribution for the bursts from SGR J0501+4516.

Zoom In Zoom Out Reset image size

We compare the two duration distributions (in count and photon space) in Section 3.4.

The emission time τ₉₀ (τ₅₀) for each burst was determined by adding the time bins (2 ms each) of high fluence in decreasing fluence rank until 90% (50%) of the fluence was reached. The emission time interval, thus, characterizes the duration of high-fluence emission. This parameter was also first introduced for GRBs as complementary to their T₉₀ (T₅₀) duration measures (Mitrofanov et al. 1999). Panels (c) and (d) in Figure 2 show the distributions of τ₉₀ and τ₅₀. These are also fit with log-normal functions obtaining 〈τ₉₀〉 = 70.3^+7.2_{− 6.5} ms (σ = 0.39 ± 0.04) and 〈τ₅₀〉 = 20.9^+2.5_{− 2.3} ms (σ = 0.30 ± 0.05).

The ratio δ₉₀ = τ₉₀/T₉₀ (δ₅₀ = τ₅₀/T₅₀) is defined as the duty cycle of each burst by Mitrofanov et al. (1999). This value should not be over 1, because the emission time excludes low fluence intervals in the burst, which sometimes are included in the duration. Panels (e) and (f) in Figure 2 present the distributions of δ₉₀ and δ₅₀, which were fit with normal distributions with 〈δ₉₀〉 = 0.68 ± 0.03 (σ = 0.14 ± 0.03) and 〈δ₅₀〉 = 0.68 ± 0.02 (σ = 0.19 ± 0.02). Both duty cycle distributions have the same mean value of 0.68.

The durations in photon space are estimated using the deconvolved source spectra with the response of the GBM detectors. In principle, since the spectral energies where SGR bursts emit most of their photons are relatively narrow (∼8–100 keV), there should not be a large difference between photon and count durations. Figure 5 shows T^ph₉₀(T₅₀^ph) versus T₉₀(T₅₀). We note that in general the count space T₉₀ values tend to be larger by a very small amount, mainly due to the fact that they take into account number of counts irrespective of their energy content. Since SGRs have mostly soft spectra, the counts corresponding to the lowest energy photons often do not contribute much flux in the durations (i.e., less than 5%). The T₅₀ estimates, however, are perfectly aligned along the x = y line, making these spectrally independent measures for durations (as also noted for GRBs by Bissaldi et al. 2011). It is also important to note here that these photon durations validate earlier duration measurements, which have been done in count space for all magnetar candidates.

Zoom In Zoom Out Reset image size

We now compare in Figure 6 the mean values of the count space T₉₀ and τ₉₀ of SGR J0501+4516 with those of four other magnetar candidates. SGR J1550–4518 was also observed with GBM and the temporal parameters are estimated using exactly the same procedures as here (van der Horst et al. 2011, in preparation). SGRs 1806–20 and 1900 + 14 were estimated using RXTE/Proportional Counter Array (PCA) observations in 2–60 keV by Göğüş et al. (2001). AXP 1E2259 + 586 durations are from Gavriil et al. (2004) and are also estimated using RXTE/PCA observations in 2–60 keV. It is obvious from the figure that all durations fall well within the same order of magnitude (∼100–150 ms), indicating a similar origin for the bursts across the magnetar population.

Zoom In Zoom Out Reset image size

Using the pulse ephemeris reported in Göğüş et al. (2010), we have aligned each burst peak to the pulse phase of the spin period of SGR J0501+4516, to search for possible correlations of the burst activity with rotational phase. In an effort to confirm the GBM barycentric correction for this analysis, we compared the barycenter corrected light curves of burst bn080823.020 (Table 1) for the GBM and RXTE/PCA data. Note that the pulse ephemeris is based primarily upon RXTE/PCA data. A cross-correlation of the GBM and RXTE/PCA burst time histories indicates no significant shift with an upper limit of 4 ms. Given the 5.76 s pulse period of this SGR, even a 4 ms shift is negligible. As a final check, the barycentering software was also tested by epoch folding TTE data to obtain a pulse profile of the Crab pulsar using the Jodrell Bank ephemeris (http://www.jb.man.ac.uk/pulsar/crab.html). We found that the phase of the first peak agreed within 200 μs of the RXTE result shown by Rots et al. (2004).

Burst peak phases were computed using the event times recorded in the TTE data in the energy range 8–60 keV. For each event time history, the burst peak time was defined as the average event time for the six most closely spaced counts in a 4 s interval surrounding the trigger time. The phase of pulse maximum was defined by fitting an inverted parabola to the folded RXTE/PCA 2–10 keV light curve. Figure 7 shows the distribution of the phase offsets of all burst peaks relative to the pulse maximum. The average offset (average of the absolute value) is 0.285 cycles. For a random phase distribution (i.e., null hypothesis), one would expect an average offset of 0.25 cycles. A Monte Carlo simulation of 10⁵ realizations for 29 draws from a random distribution shows that the probability of getting an offset value of 0.285 for 29 samples is ∼26% or roughly 1σ. We conclude that there is no correlation of the SGR J0501+4516 burst peaks with pulse phase.

Zoom In Zoom Out Reset image size

We fitted the time-integrated spectra of all 29 bursts with several models: a single power law (SPL), an optically thin thermal bremsstrahlung (OTTB), a single blackbody (BB), a power law with an exponential cutoff (COMPT)²⁰, a two blackbody spectrum (BB+BB), a single blackbody with a power law (BB+PL), a single blackbody with an OTTB (BB+OTTB), and finally, a single blackbody with an exponential cutoff power law (BB+COMPT). The parameters of the last two models (BB+OTTB, BB+COMPT) could not be constrained by most burst data, while an SPL was always a bad fit; all three models were, therefore, rejected from further spectral analysis. To determine the goodness of fit for the remaining models we used the Castor-modified Cash-statistic (C-stat). This is a modified maximum likelihood estimator which asymptotes to χ², used when there are small numbers of counts/bin (Poisson regime), which is the case for most of the SGR events (especially in the higher energy bins).

Using C-stat we were able to further reject single BB, OTTB, and BB+PL models. The first two models fit only part of the weaker burst set; for these the COMPT model also did not give a significantly smaller C-stat. It was not excluded, however, because this model, contrary to the first two, fit all bursts. Moreover, a COMPT index of 1 or −1 reproduces the BB or OTTB spectral shape. The BB+PL model overall had worse C-stat values compared to the remaining two models (COMPT and BB+BB). The relative goodness of fit among these models is exhibited in Figure 8, which shows spectral fits with COMPT, BB+BB, OTTB, and a single BB of one bright burst from SGR J0501+4516 (bn080826.136). From the figure, we see that the COMPT and BB+BB models can fit the data equally well; the residuals, however, of the OTTB and the BB model fits are unacceptably large.

Zoom In Zoom Out Reset image size

The COMPT model fits all 29 bursts well, with the BB+BB model giving equally good fits in only 18 events, where we have enough statistics to constrain the model parameters. To determine whether the COMPT or the BB+BB model fit the data best, we simulated (using RMFIT) a large set of bursts with different intensity and spectral shape parameters, using COMPT and BB+BB as input models, and then fit them with both the COMPT and BB+BB models (see van der Horst et al. 2011, in preparation, for a detailed description of the simulations). We show that the C-stat improvement is not significant enough to conclude that BB+BB (with one more parameter) is better than COMPT. A set of simulations using Xspec gave similar results to RMFIT. Below we discuss our COMPT model fits for all 29 bursts. Since the BB+BB fits provide significant information on the source parameters (Olive et al. 2004; Israel et al. 2008), we also describe the results of these fits for 18 bursts. Section 5 expands on the importance of each model.

4.1.1. COMPT Model Fits

All SGR J0501+4516 bursts were well fit with the COMPT model. We list the model parameters and statistics in Table 1 (Columns 7–11). The distributions of index and E_peak are displayed in Figure 9. The top panel shows the spectral index distribution, which is centered around zero and is best fit with a normal function with an average of −0.32 ± 0.11 (σ = 0.9 ± 0.1). This index distribution clearly excludes a pure OTTB or BB fit for all bursts, as such fits would require indices of ∼ − 1, ∼ + 1, respectively. The bottom panel of Figure 9 exhibits the E_peak distribution, which also follows a normal function with mean at 39.8 ± 0.9 keV (σ = 9.0 ± 1.0 keV). In both panels, the hatched areas highlight the distribution of 18 bursts also fit with the BB+BB model. This subsample was fit with a normal distribution with 〈E_peak〉 = 36.5 ± 1.6 keV, σ = 6.1 ± 1.4 keV and 〈index〉 = −0.63 ± 0.04, σ = 0.62 ± 0.04. Most bursts here have negative photon index and lower E_peak, indicating a softer spectrum.

Zoom In Zoom Out Reset image size

Using the COMPT spectral fits we estimated the event fluences (8–200 keV, also listed in Table 1) and plotted in Figure 10 their correlation with spectral indices (top panel) and with E_peak values (bottom panel). We note that the brightest events have a constant index value of ∼0, which progresses to lower (softer) values for weaker ones and is widely scattered with larger errors for the faintest events. Similarly, a simple trend cannot describe the relation between the hardness (E_peak) and fluence. For bursts with high energy fluence (e.g., >5 × 10⁻⁷ erg cm⁻²), $\textit {E}_{\rm {peak}}$ values are constant and cluster around 35 keV, while the values of weaker bursts cluster around 45 keV with a large scatter range between 30 and 60 keV. We note here the large span of fluence in our data (3 × 10⁻⁵ to 3 × 10⁻⁸ erg cm⁻²), which is one order of magnitude broader than earlier results (Fenimore et al. 1994; Göğüş et al. 2001; Gavriil et al. 2004). We expand on these trends in Section 5.

Zoom In Zoom Out Reset image size

We now proceed to estimate the cumulative energy fluence during the active period of SGR J0501+4516, shown at the top panel of Figure 11. The plot levels off at 3.4 × 10⁻⁵ erg cm⁻², which should be considered a lower limit (as we have not taken into account event saturation, untriggered events without TTE data, and additional bursts seen, e.g., with Swift). Assuming that the source is 2 kpc away from the Earth, then this fluence would correspond to a total energy of at least 1.8 × 10⁴⁰ erg emitted from the magnetar in bursts during this 13 day active period. The bottom panel in Figure 11 presents the differential distribution of energy fluence (log N −log S). The best fit with a power-law function is overplotted in the figure. The index of the power law is −0.48 ± 0.02, which corresponds to dN/dF∝F^−1.48. This slope is similar to the one estimated with the RXTE/PCA for SGR 1806–20 (Göğüş et al. 2000), but differs from all other SGR and AXP slope estimates, which are all very close to −1.7 (Woods & Thompson 2006). It also differs from the slope estimate of the SGR 1806–20 events detected with the International Cometary Explorer and the Burst And Transient Source Experiment data, which is −1.7 as well (Göğüş et al. 2000). However, the GBM fluence range covers the higher end of fluence values, and the analyzed burst sample is most likely incomplete at the lower fluence end (which contributes most of the events in all other sources).

Zoom In Zoom Out Reset image size

4.1.2. BB+BB Model Fits

Eighteen bursts were spectrally fit equally well with the BB+BB model. We also chose to present these results here, although this model is more complicated than the COMPT model (with one more parameter), because it has been successfully used in the past for the spectral analysis of SGR bursts (Olive et al. 2004; Feroci et al. 2004; Israel et al. 2008); we will compare our spectral results in Section 5. Assuming that the two BB components arise from two hot spots on the surface or photosphere of the source, each emission area, R², can be calculated from the temperature, T, of the corresponding BB spectrum as

$$\begin{equation} R^{2}=FD^{2}/\sigma T^{4}, \end{equation} \tag{ 3 }$$

where F is the average energy flux per event (total burst fluence divided by the spectral integration time), D is the distance to the magnetar (assumed to be 2 kpc), and σ is the Stefan–Boltzmann constant.

We plot the emission area of both BB components as a function of kT in Figure 12. As shown earlier by Israel et al. (2008) for SGR1900 + 14, we also note a clear separation of the two temperatures and emission areas of the cooler and hotter blackbodies. The cooler BB has a larger emission area, and the temperature spreads between 3 and 7 keV. The emission area of the hotter BB with temperatures in the 10–20 keV range is much smaller. The emission areas of both BB components have similar evolution through time during the active period of SGR J0501+4516 (Figure 13). The relative increase in emission area, however, is higher (by a factor of 10) for the hotter BB component. Finally, we note in Figure 14 that the total fluence of the burst is divided equally between the two BB components (we find a correlation coefficient of 0.94 corresponding to a chance probability of 5.21 × 10⁻⁹). An SPL with an index of 1.00 ± 0.05 fits the data well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We performed time-resolved spectral analysis with 8 ms temporal resolution for the five brightest bursts of SGR J0501+4516. We binned the data requiring a significance of at least 3σ above background for each bin. Each bin was then fit with the COMPT model, which was the best model for the time-integrated spectra. Figure 15 displays the light curves of these five events overplotted with their E_peak values (left column) and the correlation between E_peak and energy flux (right column). We see that the E_peak follows the light curve for the brightest part of the burst in four out of five cases (except for the fourth panel from the top), and it rises surprisingly at the beginning and the tail end of each event.

Zoom In Zoom Out Reset image size

To investigate this trend, we plotted in Figure 16 the combined E_peak values from all five bursts versus their fluxes (top panel). We clearly see that E_peak rises at both high and low flux values with a minimum determined with a broken power-law fit at 8.7 ± 0.9 × 10⁻⁶ erg cm⁻² s⁻¹. The same trend is seen in the lower panel between E_peak and fluence with a minimum at 3.4 ± 0.2 × 10⁻⁷ erg cm⁻². The indices below and above this critical flux (fluence) values are −0.28 ± 0.03 and 0.14 ± 0.01 (−0.14 ± 0.01 and 0.28 ± 0.02), respectively. An F-test comparison of a broken power law and an SPL fit shows that the former provides a better description of the correlation, with a small probability of chance coincidence 1.2 × 10⁻⁷ for the flux and 9.8 × 10⁻⁷ for the fluence. We also performed a Spearman rank test for both branches of each plot and found a low flux (fluence) correlation coefficient of −0.66 (−0.59) with a probability of chance coincidence of 1.7 × 10⁻⁹ (2.7 × 10⁻¹¹). The correlation coefficient and probability for the high flux (fluence) parts are 0.51 (0.78) and 8.7 × 10⁻⁶ (3.3 × 10⁻⁶), respectively. The latter probabilities are less significant than the ones for the low flux/fluence parts, because the Spearman rank tests involve a smaller number of data points. For the full flux (fluence) data sets, the correlations are not very significant, with a coefficient of −0.23 (−0.28) and probability of 6.4 × 10⁻³ (1.3 × 10⁻³), giving further support to changing trends (minima) in the correlations between E_peak and flux/fluence.

Zoom In Zoom Out Reset image size

Figure 17 displays the distribution of the COMPT index and E_peak of all time-resolved spectra. Both distributions were fit with normal functions with 〈E_peak〉 = 33.7 ± 0.5 keV (σ = 8.5 ± 0.5 keV) and 〈index〉 = −0.11 ± 0.09 (σ = 0.47 ± 0.09). These values are very similar to those obtained with the time-integrated fits.

Zoom In Zoom Out Reset image size

Finally, we fit the time-resolved spectra with a BB+BB model. Retaining only the fits which constrained the BB+BB parameters, we calculated the emission areas for each BB component. Similar to the integrated spectra, Figure 18 shows that the emission areas follow different behaviors with temperature for the hard and soft BB. Interestingly, a comparison with similar results shown in Figure 5 of Israel et al. (2008) for SGR1900 + 14 reveals that although the kT values are very close for the two BBs, the emission areas differ by at least one (and maybe two) order(s) of magnitude. These results indicate that either the two sources have different burst energetics or that (very unlikely) the distance determination of SGR1900 + 14 is off by a large factor.

Zoom In Zoom Out Reset image size

Using the 2 kpc distance for the source, we estimate the isotropic luminosity of the two BB components shown in Figure 19. The Spearman rank correlation coefficient is 0.83, corresponding to a chance coincidence probability of 2.78 × 10⁻¹², and indicating that these components are well correlated, as also shown with a power-law fit with index 1.1 ± 0.1. Israel et al. (2008) performed a similar analysis for the data of SGR1900 + 14 and found a spectral index of 0.70 ± 0.03.

Zoom In Zoom Out Reset image size

Our spectroscopic fitting clearly indicates that the COMPT and BB+BB models yield superior fits to the other possibilities. The COMPT model, with its power-law shape curtailed by an exponential turnover, is intended to mock-up the classic unsaturated Comptonization spectrum realized in models of accretion disks such as in Cyg X-1 or in active galactic nuclei (see Chapter 7 of Rybicki & Lightman 1979 for a summary of its development as a solution of the Kompaneets equation). These models use hot, thermal electrons in a corona to repeatedly scatter low-energy photons, heating them gradually up to an energy E ∼ kT_e consistent with the electron temperature T_e, at which point a quasi-exponential spectral turnover emerges as further heating becomes impossible. The power law marks the scale-independence of the Compton upscattering, and its slope depends only on the mean energy gain per collision (〈ΔE〉 = 4kT_e for non-relativistic electrons) and the probability of loss of photons from the scattering zone, being expressed as

$$\begin{equation} \frac{\it dN}{\it dE}\;\propto \; E^{\lambda },\quad -\lambda \; =\; -\frac{1}{2} + \sqrt{\frac{9}{4} + \frac{4}{y}}, \end{equation} \tag{ 4 }$$

with the Compton y-parameter in the domain y < 1. This parameter is the product of the average fractional energy change per scattering and the mean number of scatterings and is 4kT_e/(m_ec²) max{τ, τ²} for a non-relativistic situation for repeated Compton upscattering. Here τ is the scattering optical depth.

In the context of magnetars, such Comptonization can also ensue, but the strong magnetic field now plays an important role. Ephemeral coronae of hot electrons could be expected in a dynamic inner magnetosphere. For example, this could be due to intense dissipation of magnetic energy in the closed field line region via field line twisting, i.e., transient departures from poloidal geometry, as in the considerations of Thompson et al. (2002), Thompson & Beloborodov (2005), Beloborodov & Thompson (2007), and Nobili et al. (2008). Such a pumping of energy into electrons in low altitude regions would then be subject to irradiation by the intense bath of surface X-ray emission. The electron coronae would mimic those of their black hole counterparts discussed above and serve as a Comptonization target for the X-rays. Temporally, the coronae could be quite variable, resulting in varying or chaotic time profiles such as are observed. If the deposition of energy in hot electrons is persistent over many light-crossing timescales, the upscattered hard X-ray emission can be also. In this magnetic case, the spectrum would again naturally be a power law truncated by an exponential tail, whose energy is pinned by T_e. The slope does differ somewhat in value from the non-magnetic case because the presence of the intense field anisotropizes the environment. This alters the average fractional energy change per scattering from the isotropic 4kT_e form in Rybicki & Lightman (1979). We note that this picture is similar to that of Lyutikov & Gavriil (2006; see also Rea et al. 2008) who modeled the steep X-ray tails below 10 keV in quiescent magnetar emission using a resonant cyclotron scattering picture.

Since the electrons will move along the field lines in the zeroth Landau level, the exact kinematics that impact the determination of 〈ΔE〉 depend on the colatitude and altitude of the collisions. Therefore, the effective magnetic Compton y-parameter, y_B, that could be substituted in Equation (4), would take a somewhat different value from its non-magnetic cousin, but one anticipates that the range of spectral realizations would be similar to the B = 0 case. Since head-on collisions generally yield greater heating of photons, and these are more likely at lower altitudes, it is expected that interaction zones close to the stellar surface would spawn larger y_B and therefore flatter photon spectra. If there were a coronal radius expansion, this would then translate to an evolutionary steepening of the spectra with time. Note that the photon retention probabilities in the coronae will need to be higher than in the Lyutikov & Gavriil (2006) scenario to generate the requisite moderate optical depths τ to match the flat outburst spectra discussed in this paper.

The scatterings will take place below the cyclotron resonance unless their altitude is above around 10 stellar radii (e.g., Lyutikov & Gavriil 2006). Below resonance, the cross section is far inferior to the Thompson value (e.g., see Herold 1979). Then, the coronal electron density must accordingly be much higher than when B = 0 in order to establish a sizable optical depth. This provides a possible distinction between the two classes of magnetar emission: the steady <10 keV signal may originate at higher altitudes where the resonance is accessed and the local electron density is low, whereas the bursts reported here may be triggered closer to the surface in higher density zones that initially precipitate scattering below the cyclotron resonance. As photons get energized above 20 keV, the altitude where the resonance is accessed is lowered, increasing the cross section for scattering and therefore y_B. Note also, that to a considerable extent, polarization mode-switching (e.g., Miller 1995) between the two photon polarizations can help increase the opacity, adding a further nuance to consider when modeling the scattering environment.

Observationally, it is difficult to discern unambiguously the presence of the intense field using a Comptonization model scenario in the energy range of data presented in this paper: details of both the field strength and the geometry are subsumed in a single parameter y_B. Since the emission in this scenario should be strongly polarized, and the degree of polarization should depend on the interaction geometry, a hard X-ray polarimeter would provide insightful probes into the presence of a strongly anisotropizing super-critical field.

The dual BB can also be envisaged as a viable alternative from a theoretical standpoint. The moderate Thompson depths required to generate flat Comptonization spectra could easily be higher, thereby pushing the electron–photon interactions more toward an equilibration. The saturation temperature is then controlled by the largely uncertain total energy dissipated in the inner magnetosphere per hot electron present. Since any equilibration will be non-uniform over a coronal volume, there should be a modest temperature gradient throughout, smearing out the continuum. The range of temperatures will not be great because the system is not gravitationally hydrostatic in character. Moreover, radiative transfer effects impact the spectral shape and further modify it. Accordingly, pure BB shapes are not expected. It is quite conceivable that a two-component BB fit may well approximate the emergent continuum that is a superposition of distorted blackbodies spanning a small range of temperatures. Most probably, due to general thermodynamic considerations, the base of the coronae (e.g., centered near the source of magnetic dissipation) should be hotter than the outer layers (see also Lyubarsky 2002 and Thompson & Duncan 1995). Interestingly, the fits here generate a smaller volume associated with the hotter BB contribution, consistent with the expectation for coronal structure. Yet, somehow, we must obtain a view of the hotter zone, therefore indicating a strongly aspherical coronal geometry. Modeling this semi-equilibration is a challenging task for theorists considering the influences of the field, the twisted magnetospheric geometry, and the inherent anisotropy and polarization-dependence of the scattering process. At present, it is not possible to distinguish between this thermal scenario and a Comptonization one.

Since the COMPT model fits all events, we used it to derive fluences (time-integrated and time-resolved) and determined—using the time-resolved values—that the hardness–fluence correlation can be described by a broken power law, with a minimum at a fluence of 1.7 × 10⁻⁷ erg cm⁻². We have used E_peak values to characterize the hardness of the events in our sample. Earlier studies of this correlation used hardness ratios to bypass proper spectral fits due to low quality or insufficient data. Moreover, these studies (Fenimore et al. 1994; Göğüş et al. 2001) had a very narrow overlap in fluence space. Fenimore et al. (1994) used 95 SGR 1806–20 events detected with the Interplanetary Sun Earth Explorer-3 (ISEE-3) with fluences ranging between 1.25 × 10⁻⁷ and 8 × 10⁻⁶ erg cm⁻² and found a slightly positive correlation between the two quantities (hardness increasing with fluence). On the other hand, Göğüş et al. (2001) analyzed 159 and 385 events from SGRs 1806–20 and 1900 + 14, respectively, observed with RXTE, with fluences ranging between 1.0 × 10⁻⁹ and 2.0 × 10⁻⁷ erg cm⁻², and found the opposite trend (an anticorrelation between hardness ratios and fluences). Since the GBM data sample covers both fluence ranges (2.0 × 10⁻⁸–2.0 × 10⁻⁵ erg cm⁻²) we were able to establish that both trends are indeed correct and to define the turning point in the E_peak−fluence diagram.

Only two magnetar candidates were observed with GBM to emit a multitude of bursts, thus allowing us to construct their E_peak–flux diagrams. We find for the time-resolved data of SGR J0501+4516 that E_peak reaches a minimum of ∼30 keV at a flux value of 8.7 × 10⁻⁶ erg cm⁻² s⁻¹. The second source (SGR J1550–5418; van der Horst et al. 2011, in preparation) has flux values of 4.4 × 10⁻⁶ erg cm⁻² s⁻¹ at a similar E_peak minimum. In addition we used the hardness ratio–count rate relationship of SGR 1806–20 found in the INTEGRAL data (Götz et al. 2004) to derive an approximate flux value of 1.7 × 10⁻⁶ erg cm⁻² s⁻¹ (using their conversion from count rate to flux) at the hardness ratio minimum (Figure 3 in Götz et al. 2004). We then converted these values into isotropic source luminosity, L_iso (see Table 3); we note that (with the caveat of small number statistics and uncertainties in the distance measurements) these values are comparable [(0.4–1.5) × 10⁴⁰ erg s⁻¹], although the B-fields and fluxes at E_peak minima vary by approximately a factor of 10 between the two GBM sources and SGR 1806–20. Whether these differences reflect intrinsic source properties or instrumental effects is not yet clear. The E_peak trend is, however, clearly established in at least two different instrument data sets. The physical interpretation of this trend is beyond the scope of this paper.

Part of our sample (18 bursts) was fit with the BB+BB spectral function applied to time-integrated intervals; the remaining 11 events could not be fit due to poor statistics. In addition for five bright events we performed time-resolved spectral analysis. We found that the temperatures, emission areas and fluences (luminosities) of the two thermal components exhibit very similar properties in both cases (time-integrated and time-resolved). Our results are consistent with those presented for intermediate and short bursts from SGR 1900 + 14 by Olive et al. (2004) and Israel et al. (2008); the SGR J0501+4516 burst emission areas and temperatures fall into the same region with those of the short bursts of SGR 1900 + 14 shown in Figure 5 of Israel et al. (2008). We also see a similar trend in the correlation between the luminosities of the two components as described by Israel et al. (2008), namely that both luminosities increase in tandem. This behavior indicates that the hot and cool BB components may come from two separate emission regions, as also pointed out by Israel et al. (2008): a smaller but hotter one from the surface of the magnetar, and a larger, cooler one from the star's magnetosphere (but see also the discussion in Section 5.1). We note here that the time-integrated hot BB emission area of 13 SGR J0501+4516 bursts (i.e., not including the five brightest events; Figure 12) is similar to the emission area (∼0.05 km²) of the BB component found in the persistent emission of SGR J1550−5418 during one active bursting episode in 2009 January, believed to originate from a hot spot on the neutron star surface (Kaneko et al. 2010).

Finally, we performed a detailed temporal analysis of all 29 bursts and estimated their durations (T₉₀/T₅₀) for the first time in count and in photon space. Both estimates agree within statistics and thus validated all earlier (count space) duration estimates of magnetar candidate bursts. The durations (and for four SGRs also the emission times and duty cycles) of five magnetar candidates follow a log-normal or normal distribution. We find that SGR J0501+4516 events are very similar in average duration with four more magnetar candidates (three SGRs and one AXP). However, the T₉₀ distribution of the bursts from AXP 1E2259 + 586 has a factor-of-two larger dispersion (σ ∼ 0.73, Gavriil et al. 2004) than those of at least three other SGRs (σ ∼ 0.34, 0.35, Göğüş et al. 2001, and ∼0.35, present work). SGR J0501+4516 bursts have an average duty cycle (δ₉₀ = 0.68) larger than other SGRs (0.45, 0.46 for SGRs 1900 + 14, 1806–20, respectively; Göğüş et al. 2001). The differences in the intrinsic properties of the sources may be due to differences in the size of the active region responsible for the burst emission. We discuss below the burst energetics and its evolution. Similarly to other magnetars, we also do not find a correlation between the pulse phase of the source and the burst peak times.

The overall behavior of SGR J0501+4516 during its 13 day active period (2008 August 22 to September 3) is also very interesting. Although about half of the bursts (16/29) occurred on one day (2008 August 23), most of the burst energy was emitted from the source during August 24–26 (see also Figure 11). The average fluence of the beginning and the later part of the active period was constant and at a lower level. During these three days, GBM detected seven bursts, five of which are also distinguished in that the emission areas from both BB components are the largest. Further, one of these five events (bn080824.054) shows a double-peaked structure reminiscent of some bright thermonuclear Type I X-ray bursts from accreting neutron stars, which has been interpreted as due to photospheric radius expansion (PRE). Watts et al. (2010) studied the effects of PRE in high magnetic fields using this event and find that the predicted flux from PRE theory is consistent with the one observed, opening the way to determining fundamental parameters of neutron stars such as their equation of state.