QUASI-PERIODIC OSCILLATIONS IN SHORT RECURRING BURSTS OF THE SOFT GAMMA REPEATER J1550−5418

Soft gamma repeaters (SGRs) represent a small class of neutron stars whose slow spin periods and high spin-down rates imply an unusually strong magnetic field in the excess of 10¹⁴ G. Duncan & Thompson (1992) and Thompson & Duncan (1995) predicted the existence of such objects, which they named magnetars. SGRs are believed to be one of two observational manifestations of neutron stars with an exceptionally strong magnetic field; anomalous X-ray pulsars (AXPs) form the other class of objects, although evidence suggests that there is no clear-cut line between them, and recently a low magnetic field source has been found (Rea et al. 2010).

The defining characteristic of SGRs are irregular bursts that vary in duration from tens to hundreds of milliseconds and span ∼5 orders of magnitude in peak luminosity (10³⁸ to 10⁴³ erg s⁻¹) in hard X-rays <100 keV. However, there is a very rare type of burst, the so-called giant flares, which have been only detected three times in the past 34 yr from three different sources. These reach peak luminosities of ∼10⁴⁵ erg s⁻¹ and are believed to be powered by a catastrophic reordering of the magnetic field. Since this field is coupled to the solid crust, Duncan (1998) suggested that such large-scale reconfiguration might rupture the crust, triggering global seismic vibrations that would be visible as periodic modulations of the X-ray and γ-ray flux. This idea was confirmed by the detection of quasi-periodic oscillations (QPOs, i.e., stochastic processes that vary on a characteristic timescale) in the expected range of frequencies (∼10–1000 Hz) in the tails of giant flares from two different magnetars (Israel et al. 2005; Strohmayer & Watts 2005, 2006; Watts & Strohmayer 2006). SGR giant flares thus present outstanding test cases for testing theories of neutron star structure and composition models. Several intermediate flares, in energy and duration between the short bursts and the giant flares, have also been observed, but no QPOs have been found in these bursts (Watts 2012).

To date, there have been few searches for QPOs in recurrent bursts of magnetars. El-Mezeini & Ibrahim (2010) reported QPO detections in a sample of bursts from SGR 1806−20 observed between 2 keV and 60 keV with the RXTE, however, a revised analysis incorporating variability in the burst envelope showed that the reported QPOs are not significant (Huppenkothen et al. 2013).

Finding QPOs in short SGR bursts is technically challenging: as shown in Huppenkothen et al. (2013), standard Fourier methods commonly used for this purpose fail when applied to the short, highly variable burst light curves. The major difficulty lies in the non-stationarity of magnetar bursts. The statistical distributions generally used in Fourier analysis in astronomy are strictly only valid for processes whose properties do not vary over the duration of an observation. This is clearly not true for an SGR burst: they are short events, exhibiting variability on timescales roughly equivalent to the periods of QPOs observed in the giant flares. Below 100 Hz or so, many of the bursts exhibit a wealth of variability properties: to leading order, there is the rise and fall of the burst, i.e., a burst envelope. In most bursts, the envelope has a high degree of temporal substructure beyond this envelope. This substructure differs widely from burst to burst and is poorly understood. Perhaps what we call a burst is actually a superposition of many smaller events. Alternatively, the overall burst shape could be composed of an envelope combined with a stochastic process, leading to additional variability on shorter timescales. Finally, the complexity of the burst envelopes in general varies with energy, inserting another constraint in our interpretation of their structure. This lack of knowledge leads to two major problems when searching for QPOs. At low frequencies, very few cycles of a potential oscillation are sampled due to the short duration of the burst. A succession of peaks may look like a quasi-periodic signal to the naked eye, but could be a chance superposition of a stochastic process, without the characteristic timescale implied by a QPO. The other major difficulty is our lack of knowledge of the underlying statistical distribution that we must test against. The statistical distributions generally used in testing for QPOs are strictly defined for stationary stochastic processes. The presence of a burst envelope changes the observed distributions at low frequencies from those we know. This makes it difficult to derive inferences about the presence of a QPO at these frequencies.

In the absence of this knowledge, it is possible to make reasonable assumptions. In Huppenkothen et al. (2013), we introduced a Bayesian approach to deal with our uncertainty in the underlying burst processes by assuming a purely stochastic process with a power-law power spectrum, a so-called red noise process. While this assumption is strictly not true, either, we showed that it is a conservative choice: in practice, the presence of the burst envelope narrows the statistical distributions at low frequencies compared to the distribution we use to model the process. We are thus more likely to underestimate the significance of a signal at low frequency than overestimating our confidence in a detection. In Huppenkothen et al. (2013), we also analyzed a short bursting episode of the magnetar SGR J0501+4516, where we found one candidate detection out of 27 bursts. Our results were inconclusive with regard to the origin of this signal and showed where our method can potentially produce ambiguous results: the significant detections were all at integer multiple frequencies of a low-frequency signal at ∼30 Hz, which was heavily affected by red noise and thus not significant itself. However, this signal corresponds to less than two full cycles at 30 Hz, given the short duration of the burst. We thus concluded that this signal was equally likely to be a chance occurrence of two red noise peaks close together as it was to be a QPO, and we deferred a more in-depth discussion to a later work with a larger sample of bursts.

In this paper, we perform a comprehensive search for QPOs in a much larger sample of bursts from a so-called burst storm observed from SGR J1550−5418 in 2009 January. SGR J1550−5418 (also 1E 1547.0−5408) was first observed with the Einstein X-Ray Observatory (Lamb & Markert 1981). Later observations with XMM-Newton revealed a soft X-ray spectrum and a possible association with a young supernova remnant, suggesting that it might be an AXP (Gelfand & Gaensler 2007). The AXP nature was confirmed by the subsequent detection of radio pulsations with a slow spin period of P = 2.096 s and a spin-down of $\dot{P} = 2.318 \times 10^{-14}$, implying a magnetic field of 3.2 × 10¹⁴ G (Camilo et al. 2007).

SGR J1550−5418 exhibited three major bursting episodes: in 2008 October, 2009 January, and 2009 March/April. The January episode was exceptional: the source showed hundreds of bursts within a single day, observed with several X-ray telescopes: the Swift Burst Alert Telescope (BAT; Israel et al. 2010; Scholz & Kaspi 2011), the Fermi Gamma-Ray Burst Monitor (GBM; Kaneko et al. 2010; von Kienlin et al. 2012; van der Horst et al. 2012), RXTE (Dib et al. 2012), and two main instruments on board the INTEGRAL spacecraft (Mereghetti et al. 2009; Savchenko et al. 2010).

Burst storms like the one observed from SGR J1550−5418 are rare and have been observed in only three other sources (SGR 1806−20, SGR 1900+14, and SGR 1627−41; Götz et al. 2006; Israel et al. 2008; Mazets et al. 1999), the first two of which have also exhibited a giant flare. During the first triggered observation recorded with Fermi/GBM on 2009 January 22, the source also showed an increase in persistent flux level up to ∼100 keV (Kaneko et al. 2010) for around 150 s of intense bursting. A subsequent search for pulsations in this plateau of hard emission revealed a signal at the period of the neutron star, but no higher-frequency QPOs. The bursting episode ended in 2009 April, and there have been no subsequent bursts recorded since. A catalogue of magnetar bursts observed with Fermi/GBM is currently in preparation (A. C. Collazzi et al., in preparation).

Here we present the first large-scale robust QPO search from the 2009 January burst storm, observed with Fermi/GBM. In Section 2 we describe the sample in some detail, and in Section 3 we give a very brief overview of the Bayesian technique used in the QPO searches. Finally, in Section 4 we present our results and interpret both the QPO searches as well as a characterization of broadband variability in the bursts in Section 5.

X-ray bursts from SGR J1550−5418 triggered Fermi/GBM a total of 55 times between 2009 January 22 and 29, with 41 triggers on January 22 alone. Each trigger records data from 30 s before the trigger up until 300 s after the trigger. As a result, multiple (untriggered) bursts were observed in most triggers. van der Horst et al. (2012) identified a total of 286 bursts in this sample, which have time-tagged event (TTE) data available. The TTE data type has a time resolution of 2 μs, needed for high-precision timing studies. We use data from the 12 Na i detectors, whose energy range of 8 keV to 4 MeV is sufficient, since SGR bursts rarely exhibit radiation above 200 keV. Additionally, we only used detectors with viewing angles to the source <60° and checked whether the source was occulted by the spacecraft and the other instrument, the Large Area Detector (LAT).

The sample is the same as in van der Horst et al. (2012), and we use the burst durations, start times, and fluences from that paper in our analysis. We extracted TTE data between 8 keV and 200 keV around each burst, starting at t_start − 0.1 × T90 (the burst duration, T90, is defined as the time in which the central 90% of the photons, starting at 5% and ending at 95%, reach the detector) and ending at t_start + 1.1 × T90 in order to ensure that the entire burst is within our data set. The sample has a mean duration of 0.174 s and an overall asymmetric shape with a faster rise than decay. The estimated fluences range from 10⁻⁸ to 10⁻⁵ erg cm⁻². For a more detailed description of the data extraction process and sample definition, see van der Horst et al. (2012). An analysis of the first trigger, including a timing analysis of the inter-burst periods, was performed in Kaneko et al. (2010). Time-resolved spectroscopy and the spectral evolution with burst flux are discussed in Younes et al. (2014). Of the bursts in the sample, 23 have saturated parts, where the detector cannot record all photons during periods of very high count rates. We excluded all 23 bursts from our analysis due to the rather complicated effect saturation has on the timing analysis. This gives us a total sample of 263 bursts.

QPOs are generally found by taking the Fourier transform of a light curve and looking for variability focused at a particular frequency. The square of the Fourier transform of the data is called the periodogram.¹⁰ Different types of variability have different frequency distributions. Our task becomes to disentangle the different components in the periodogram. While pure photon counting noise has a flat power spectrum with a well-behaved χ² distribution with two degrees of freedom about a constant mean, QPOs produce sharp coherent features. Stochastic processes with correlated frequencies, often termed "red noise" or "1/f noise," are also often observed and follow power laws or broken power laws with stronger variability at low frequencies and decreasing power at higher frequencies.

The short duration of magnetar bursts means that this low-frequency variability has timescales similar to those of the QPOs observed in the giant flares. Thus, we must test for QPOs in a periodogram consisting of complicated variability. We adopt the method from Huppenkothen et al. (2013), first suggested for red-noise-dominated periodograms in Vaughan (2010). This method assumes an exponential distribution of powers about the underlying power spectrum, which we assume to be a power law or broken power law. To find QPOs, we fit a broadband noise power spectrum to each burst, which is then divided out. The highest outlier in the residuals is our candidate QPO detection. We then simulate fake periodograms using the broadband noise power spectrum and incorporating uncertainties in the model parameters, and we perform the same detection procedure on those simulations. We can thus compare the observed highest data/model outlier with a distribution of data/model outliers from the simulations, to infer the probability that the observed outlier is a significant QPO.

Below, we give a very brief overview of the QPO search strategy, and we refer the reader to Huppenkothen et al. (2013) for a detailed description, discussion of the method's limitations, and tests on both simulated data and a smaller sample of magnetar bursts.

In more detail, for each burst:

• 1.  
  We fit both a power law and a broken power law, i.e., the broadband noise model, to the periodogram of the burst observation. We fit the unnormalized posterior predictive distribution, consisting of a likelihood function following a χ² distribution around the broadband noise model and priors that are independent of each other and of the form p(θ) = 1/θ (scale prior) for scale parameters (e.g., broadband noise amplitudes) and flat otherwise. As a result, we obtain the so-called maximum a posteriori (MAP) as a result of the numerical optimization step. The MAP estimate is the Bayesian equivalent of the maximum likelihood. For both models, we then construct the ratio of likelihoods at the parameter values corresponding to the MAP estimate.
• 2.  
  We sample the posterior predictive distribution of the simpler broadband noise model—the power law—using a Markov chain Monte Carlo (MCMC) technique, in this case employing an affine-invariant MCMC ensemble sampler (Goodman & Weare 2010), as implemented in python by emcee (Foreman-Mackey et al. 2013). The resulting ensemble of parameter values will follow the posterior distribution of the assumed broadband model, thus allowing for statistical inferences over this distribution.
• 3.  
  We simulate N_s artificial periodograms from the MCMC sample and fit each with the two broadband noise models considered such that we can construct the likelihood ratio for each of the fake periodograms. This will allow us to construct a distribution of likelihood ratios from a sample we know to be derived from the simpler model. If the likelihood ratio obtained for the observed periodogram is an outlier of the distribution of likelihood ratios, then the observed data are unlikely to be generated from this model. Note that this is strictly evidence against the power-law model; it is not direct evidence in favor of the more complex model. We use this approach to reject the simple power-law model for cases where the posterior predictive p-value (the ratio of samples in the posterior distribution of likelihood ratios lying above the observed values, divided by the total number of samples in this distribution) falls below 0.05. If this is true, we use the broken power law to model the broadband component of the periodogram; otherwise, the simple power-law model is adopted. A threshold of p < 0.05 is not very stringent, but desirable. We would rather reject the simpler broadband noise model in favor of a more complex one. It is preferable to overfit the periodogram, rather than underfit, because broadband noise features not adequately modeled by the broadband noise model may instead be mistaken for QPOs in the subsequent analysis.
• 4.  
  We construct a second MCMC sample from the adopted broadband noise model in the same fashion as in step 2. We simulate N_s periodograms from this sample and fit with the adopted broadband model. For the observed periodogram and each fake periodogram, we divide the periodogram by the best-fit MAP parameter estimate and define the test statistic T_R = max_j(2I_j/S_j), where I_j is the observed power I at frequency j and S_j is the value of model S at frequency j. This test statistic is the maximum power in the residual periodogram after dividing out the broadband noise model. In the ideal case where the parameters θ defining the model S are perfectly known, 2I_j/S_j follows a $\chi ^2_2$ distribution. In reality, there is an uncertainty in S_j, since the parameters θ are not known exactly, leading to a deviation in the distribution of T_R from the theoretical expectation. Sampling the posterior probability distribution of the parameters given the data and the model via MCMC allows us to construct the actual distribution of T_R from simulated periodograms, taking into account all relevant uncertainties. We can thus construct a posterior distribution for the T_R statistic under the null hypothesis that the observed maximum power is due to a stochastic aperiodic process. Comparing T_R from the real data to this simulated distribution allows us to define a posterior predictive p-value for this null hypothesis. If the latter is small, then the observed maximum power is unlikely to be a product of a $\chi ^2_2$ distribution.

Although the giant flare QPOs were all very narrow, we cannot exclude broader signals in the shorter bursts. We search for these signals by performing exactly the same analysis as in step 4 on periodograms that are binned in frequency, where the bin widths are chosen from a logarithmic scale between 1 Hz and 200 Hz. Binning a broad QPO signal makes it easier to detect, since the QPO is grouped into a single bin, while random fluctuations from one frequency bin to the next are suppressed. We search for QPOs in the same way as in the unbinned periodogram: by defining T_R for the binned periodogram and comparing to the distribution of binned T_R values from the sample of simulated periodograms. In this case the p-values for different bins are not independent. To avoid excessive false positives, we accept significant detections only if they are detected in at least two different bin widths at the same frequency. In order to conserve computation time, we set the number of simulations N_s = 10⁴. This implies that the significance can only be quoted to p = 10⁻⁴ for a single trial. The detection limit we use depends on the number of trials: the more periodograms we search, the more likely it becomes to make a significant detection purely by chance, even if no signal is present. We thus require a more stringent detection limit for searching individual bursts, where we search hundreds of periodograms, than for searching averaged periodograms, where we only search 10. For searching individual bursts as in Section 4.1, we require p < 10⁻⁴ for a single trial, corresponding to p < 0.0263 or roughly 2.3σ, given the number of bursts in our sample. For the 10 averaged periodograms we search in Sections 4.2 and 4.2.2, we choose 6 × 10⁻³, or roughly 3.5σ. All p-values given below are trial-corrected: in the search of individual bursts, we correct by the total number of those bursts, i.e., 263, and for the averaged periodograms by the number of averaged periodograms searched, i.e., 10. The number of frequencies and bin widths we searched over is automatically taken into account by our methodology, by searching over the entire frequency and bin width range for the simulations as for the real data.

We searched light curves from 263 individual bursts for periodic signals and QPOs. In order to be sure to include the entire burst, we added 10% on either side of the burst duration (T90). Additionally, we constructed averaged periodograms from samples of bursts to explore whether a signal could be re-excited in several bursts. Finally, we characterized broadband variability for the sample as a whole, which may guide future work on emission and trigger mechanisms.

4.1. Individual Burst Searches

We searched all 263 bursts for QPOs over the complete range of available frequencies from ⩽10 Hz to 4000 Hz. The maximum frequency was chosen to maximize computational efficiencies, while at the same time oscillatory modes are unlikely to occur at a much higher frequency.

Four bursts show detections significant with p < 10⁻⁴ (single-trial) or p < 0.0263 corrected for N_b = 263 trials. Three candidates are significantly affected by dead time and pileup, that is, their count rate is close to the saturation count rate. This is the case when a significant number of photons arrive within less than 2.6 μs of each other (the dead time of the GBM recoding system) and are consequently recorded as a single photon. Here we used the highest intrinsic time resolution data from GBM: TTE data with 2 μs resolution. While 2.6 μs correspond to a higher frequency than we are interested in, the above effects can nevertheless influence the periodogram in nontrivial detector-dependent ways, which are not retrievable or quantifiable. A proper treatment of affected bursts is beyond the scope of this work; we thus consider the QPO search on these bursts as inconclusive and make no further statements about their properties.

The remaining burst, one of several in TTE data of trigger 090122218 with a burst duration of 0.49 s, has a significant detection of a broad feature at 260 Hz with p < 0.0263 (trial-corrected; also includes an uncertainty in the parameters of the broadband model). We plot the light curve and periodogram of this candidate in Figure 1. While there might be some red noise power left at these frequencies, the signal is largely dominated by white noise. We use the traditional (analytical) test against white noise for an upper limit on the detection probability (Groth 1975; van der Klis 1989). This would be a precise estimate if there was no red noise in the signal, but as we cannot exclude some contamination from red noise, this must be regarded as an upper limit instead. We find that the probability that the observed peak in the periodogram is due to Poisson counting noise alone is p = 5.26 × 10⁻⁶. The fractional rms amplitude is high, rms_frac = 21% ± 3%, as estimated from integrating over the noise-level-subtracted periodogram. We estimated the error following Heil et al. (2012). This error calculation is somewhat too simplistic for the periodogram we consider here: there may be a residual contribution of aperiodic variability contaminating the powers we integrate over, which is not taken into account properly. However, our lack of knowledge about the burst processes involved precludes us from running simulations to establish the error to a higher degree of precision.

Zoom In Zoom Out Reset image size

We measured the Q-value, defined as the centroid frequency divided by the width of the signal, by comparing the p-values for periodograms of this burst binned at several frequency resolutions and picking the frequency resolution that yielded the lowest p-value to reflect the most likely width of the signal. The QPO is extremely broad: the Q-value is Q = ν₀/Δν = 2.9. This is at the lower boundary of what one would call a QPO (Q > 2) as opposed to a broadband noise feature. Due to the single occurrence of this signal in the sample, it is not possible to average periodograms, as is usually done to improve signal-to-noise ratio and estimate errors, such that a Lorentzian fit to the feature is possible for a precise estimate of the width and the rms amplitude (van der Klis 2006).

The periodogram of the same burst also shows a broad feature at 20 Hz. In order to understand the origin of this feature and its connection to the QPO at 260 Hz, we fit two fast-rise, slow-decay profiles to the light curve. We used the skew-normal distribution, a generalization of a simple Gaussian profile that allows for skewness and has the form

$$\begin{equation} f(t) = \frac{2}{\sigma } \phi \left(\frac{t-\mu }{\sigma } \right) \Phi \left(\alpha \left(\frac{t - \mu }{\sigma } \right) \right) \;, \end{equation} \tag{ 1 }$$

where

$$\begin{equation*} \phi \left(\frac{t - \mu }{\sigma } \right) = \frac{1}{\sigma \sqrt{2\pi }} \exp {-\frac{(t-\mu)^{2}}{2\sigma ^{2}}} \end{equation*}$$

and

$$\begin{equation} \Phi \left(\alpha \left(\frac{t - \mu }{\sigma } \right) \right) = 0.5 \left[ 1 + \mathrm{erf}\left(\alpha \frac{t - \mu }{\sqrt{2}\sigma }\right) \right] \;. \end{equation} \tag{ 2 }$$

Here μ is the location in time of the peak of the profile, σ is the width, and α is a skewness parameter (Azzalini 1985). We find that the signal at 20 Hz is easily reproduced by a superposition of two skewed peaks with a separation of 0.04 s and widths σ₁ = 0.02 s and σ₂ = 0.016 s. While the feature is easily reproduced by two non-periodic functions, there are too few cycles observed to make a strong statement about its nature (see Huppenkothen et al. 2013, for a similar feature). However, it cannot explain the highly significant signal at 260 Hz: the power spectrum of the two skew-normal functions fitted to the data turns over at lower frequencies and becomes negligible above 200 Hz. Beyond this frequency, there is very little power in this model, and the power spectrum at higher frequencies should be dominated by Poisson noise only. This implies that the QPO is not easily reproduced by a burst envelope and is likely a separate process producing variability at these frequencies. In order to confirm this observation, we have fit the observed light curve with both standard Gaussian profiles and Lorentzian profiles. Both alternatives give results very similar to the one presented above: a near-perfect fit to the low-frequency feature and a sharp drop in power around 200 Hz.

The sensitivity limits for signal detection vary strongly from burst to burst and with frequency, especially for the low-frequency part of the periodogram, where the contamination by broadband variability is strong. Below ∼100 Hz, sensitivities range from ∼50% fractional rms amplitude at 30 Hz to ∼10% fractional rms amplitude at 100 Hz. Above ∼150 Hz, the bursts are almost all dominated by photon detector noise, and a QPO should be the only source of non-white-noise variability in this regime. Our method converges toward standard Fourier methods in this frequency range. Instrumental effects such as dead time can still be an issue; neither method is equipped to deal with these effects without a large number of dedicated simulations. Above 150 Hz, sensitivities are generally in the range of 5%–10% fractional rms amplitude.

4.2. Averaged Periodograms

To increase sensitivity, we average the periodograms of a number of bursts. This assumes that the short bursts always excite the same star quakes, which has also been seen in giant flares, where QPOs are detected to be present over many cycles.

4.2.1. Signal Grouped by Burst Duration

We sorted the bursts by duration (T90) into five groups: <50 ms, 50–100 ms, 100–250 ms, 250–500 ms, and >500 ms. To average periodograms, we picked the longest burst in each group and extracted light curves of the same duration for each burst in the sample, so that each periodogram would have the same number of frequencies. We then averaged the periodograms within a group to get the final periodogram. Since we use light curves of equal duration within each group, the shorter bursts in each group add noise into the final averaged periodogram, which reduces the QPO detection threshold somewhat. Limiting this effect is our main reason for dividing the bursts into groups, so that we can search for QPOs in the longest bursts without a strong noise component added by including the shortest bursts in the same sample.

There are no QPOs detected in the first four averaged periodograms. We report a candidate detection in the averaged periodogram of the longest bursts (T90 > 0.5 s, N_bursts = 47; see Figure 2 for the averaged periodogram). The strongest signal with p < 2.5 × 10⁻³ occurs at 10 Hz, with a width of ∼5 Hz. Note that 10 Hz corresponds to a timescale of 0.1 s, close to the peak of the distribution of burst durations. However, we cannot exclude that this feature is actually an artifact caused by an inadequate characterization of the underlying power spectrum. Another process, such as a doubly broken power law or a combination of Lorentzians as often used in broadband noise modeling of X-ray binaries, may represent the shape of the power spectrum better, but requires more intricate model selection criteria than implemented here.

Zoom In Zoom Out Reset image size

4.2.2. Averaged Periodograms per Trigger

We also search for QPOs in trigger data sets with high-resolution TTE data. Since data sets obtained with Fermi/GBM from an individual trigger are roughly 330 s, we searched those triggers with a large number of bursts (see Table 1 for an overview) in a short time span for long-lived signals. As for the duration-averaged periodograms, we extracted light curves around all bursts in those five triggers of the duration of the longest burst in that trigger data set. We then constructed the periodograms of these light curves and computed the average periodogram of the sample. The resulting periodograms do not all have the same frequency resolution; for those with a frequency resolution less than 1 Hz, we averaged neighboring frequency bins to achieve a resolution close to 1 Hz.

Table 1. Signals from the Five Averaged Triggers, and Detection Sensitivities for Different Frequencies

Trigger ID	Nbursts	Min T90	Max T90	ν0	Δν	Posterior	Simulated	Sensitivities in Fractional rms Amplitude
		(s)	(s)	(Hz)	(Hz)	p-value	p-value	40 Hz	70 Hz	100 Hz	1000 Hz
090122037	32	0.0322	1.0724	99	27	<4 × 10−4	0.107	3.6	2.4	⋅⋅⋅	1.4
090122052	28	0.0364	1.4952	127	10	<4 × 10−4	0.016	4.8	2.5	1.9	1.5
090122194	20	0.0364	1.2124	93	12	<4 × 10−4	0.013	6.5	3.8	2.7	1.8
090122218	21	0.1176	1.3496	91	10	1.2 × 10−3	0.009	5.2	3.0	2.3	1.6
090122283	30	0.0504	2.4724	61	20	8 × 10−3	⋅⋅⋅	2.9	1.7	1.2	0.9

Notes. This table summarizes the results from the averaged periodograms of five triggers. The significant detections are shown in Figure 3. The last, 090122283, had no significant detections with p < 6 × 10⁻³, single trial probability. The second column gives the number of bursts averaged together, which equals the number of bursts in the trigger, excluding those that have saturated parts. The third and fourth columns give the minimum and maximum burst durations in the sample, respectively. Columns five and six present the centroid frequency ν₀ of the observed signal and the corresponding frequency bin width Δν in which the signal is detected, respectively. We quote the detection threshold sensitivities where no detection has been made. The posterior p-value is estimated from simulations derived from the MCMC sample of the broadband model for each periodogram. The second p-value is derived from averaging random subsets of burst periodograms and extracting the highest outlier from the data-model residuals for each averaged periodogram, as described in the text.

Download table as:  ASCIITypeset image

We searched data sets from five triggers with 20 to 32 bursts per trigger and excluded long-timescale variability below 60 Hz from the range of frequencies searched. Below 60 Hz, there will be a significant contribution from the overall shape of the short bursts (as in the 20 Hz feature discussed in Section 4.1), and thus our estimates are unreliable. We search both the unbinned periodogram and periodograms binned to different frequency resolutions between 1 Hz and 200 Hz, but considered only candidate signals with Q > 7. This is necessary, because at low frequencies, a candidate signal in a frequency bin that is wider than 0.2ν₀ likely incorporates power from the part of the periodogram below 60 Hz, where we believe estimates to be unreliable. We find candidate detections in three of the triggers. The results are summarized in Table 1. The signals, at 93 Hz (trigger ID 090122194), 91 Hz (trigger ID 090122218), and 127 Hz (trigger ID 090122052), are fairly narrow, Q ≈ 7 and Q ≈ 13, respectively. Two of the signals are significant to p < 4 × 10⁻⁴ (roughly 3.7σ), and the third has a p-value of p = 1.2 × 10⁻³ (3.4σ). A fourth signal (trigger ID 090122037) is significant to p < 4 × 10⁻⁴, but fails to fulfill our criterion of Q > 7. At the same time, this signal is at a frequency of 99 Hz, close to the frequency where significant detections were made in two of the other triggers. We plot the periodograms for all three triggers in Figure 3.

Zoom In Zoom Out Reset image size

The periodogram shape may change between different bursts, largely due to the wide spread in burst duration, fluence, and burst shape. The effects this may have on the averaged periodogram are hard to quantify without a large number of dedicated simulations of the overall burst variability, which is beyond the scope of this work. In order to test whether the observed QPOs could be due to the differences in duration, fluence, and shape in the averaged samples, we constructed a large number (N_s = 10³) of averaged periodograms from randomly selected subsets of the burst sample, excluding the four triggers where candidates were observed. If the QPOs are due to effects of the varying burst properties, then these signals should appear in a large number of these simulations. We searched these periodograms in the same way we did for the averaged periodograms from individual trigger data and compared the resulting distribution of maximum powers >60 Hz from the data-model residuals to the maximum powers from the averaged periodograms of individual triggers. Column 8 in Table 1 shows the p-values of observing the candidate signals presented above from a random subset of bursts.

In Table 1, we also show the detection sensitivities for all five averaged periodograms for various frequencies, for the ∼1 Hz resolution of the periodogram. Note that we are even more sensitive to broader signals, as averaging neighboring frequency bins increases the signal-to-noise ratio. The numbers quoted here are upper limits for the most coherent signal we could have observed and were calculated for each quoted frequency from the distribution of maximum powers of the MCMC-derived sample of periodograms. The frequency dependence of the sensitivity is due to the fact that the low-frequency part of the periodogram is dominated by aperiodic red noise variability, and any quasi-periodic signal needs to introduce strong variations in order to be visible. Sensitivities for fractional rms amplitudes are between 3% and 6% for the lowest frequencies and drop to 0.9%–1.7% at high frequencies. The differences in sensitivities between triggers are due to a combination of number of averaged bursts, number of averaged frequencies, and the average count rates of bursts included.

4.3. Broadband Variability

Magnetar bursts are a class of events with very diverse light curves: they differ vastly in duration and peak count rates, but also in overall burst shape (see Huppenkothen et al. 2013, for examples of burst light curves). How exactly this variability is produced is not well understood. Are all bursts a realization of fundamentally the same process? Are there characteristic rise or decay timescales? It is useful to characterize the variability properties of a large sample, which may answer some of these questions.

In the following, we give an overview of the broadband variability observed in the whole sample of bursts. Out of 263 burst periodograms, 193 were adequately fit with a simple power law plus a constant to account for the white-noise component; the remaining 70 rejected the null hypothesis to p < 0.05, and we thus adopted a broken power law for these periodograms. Burst duration and burst fluence could influence whether a simple power law or a broken power law fits the broadband noise. For example, for dim bursts, variability observed in the bright bursts may be hidden in the noise. In order to test this idea, we plot the fluence and burst duration distributions for bursts modeled with a power law and a broken power law in Figure 4.

Zoom In Zoom Out Reset image size

Burst duration (T90) shows only a marginal difference in the T90 distribution (p = 4.6 × 10⁻⁴ for a two-sample Kolmogorov–Smirnov (K-S) test). There is an appreciable difference in fluence between the samples (two-sample K-S test: p = 9.32 × 10⁻¹¹): bursts with broken power-law power spectra have higher fluence than bursts modeled with a simple power law. Note that the threshold for rejecting the power-law broadband model is not very high, p < 0.05. This is desirable for the main objective of our analysis, the search for QPOs: if the broadband noise is not adequately represented by the model, then broadband features may be attributed to QPOs instead, leading to false-positive detections. Setting a threshold of 5% is a compromise between reliability of QPO detection at the expense ofpotentially contaminating our sample of bursts fit with a broken power law with bursts that are consistent with a simple power law. In practice, this decreases the probability for rejecting the null hypothesis when performing the K-S test: the difference between the two distributions may be stronger than we report here.

As well as studying the overall properties of bursts with different broadband noise models, we can also study the broadband noise properties of the sample and see whether the noise properties change with fluence or burst duration. In Figure 5, we show the distribution of power-law indices for the various components. For the bursts modeled by a simple power law, the distribution of power-law indices varies between 1.5 and 4, with a median at μ_γ = 2.42. The average low-frequency component of the two-component broken power law is flatter than for a single power law ($\mu _{\gamma _0} = 1.49$), while the higher frequency component is much steeper ($\mu _{\gamma _1} = 6.16$). Note that for several bursts, the second component is extremely steep. This may be caused by the contamination of this sample with bursts that were incorrectly classified as too complex for the simple power law. In this case, the second power-law index is often not well constrained and tends to high values. In Figure 5, we also plot the break frequency between the two components of the broken power law for those burst periodograms for which the simple power law was rejected. The distribution peaks around 100 Hz, below which the power-law index is fairly flat. At higher frequencies, it steepens considerably, as shown in panel 3 of Figure 5. The distribution is fairly broad, with the bulk of burst periodograms breaking between 30 Hz and 400 Hz.

Zoom In Zoom Out Reset image size

We correlated the power-law indices with various burst properties to see whether there is a systematic effect due to burst duration or brightness, similarly to the reasoning for why some bursts require a more complex model than a simple power law. There is an indication of a correlation of the burst duration with the power-law index γ for bursts modeled with a simple power law (see Figure 6): shorter bursts seem to have slightly steeper power-law indices. A Spearman rank coefficient yields R = −0.51, indicating an anti-correlation with a probability for the null hypothesis (no correlation) of p = 8.65 × 10⁻¹⁶. In order to compute the slope of the correlation and incorporate the errors on the measurements of the power-law index, we binned the data set logarithmically into seven bins and computed the mean power-law index within each bin (following work in, e.g., X-ray binaries and active galactic nuclei, Gleissner et al. 2004). We estimated the error on the mean as a standard error, $\sigma /\sqrt{n}$, where σ is the standard deviation of the sample and n is the number of data points in each bin. The correlation can be fit with a power law, with index α = −0.1027 ± 0.00523.

Zoom In Zoom Out Reset image size

In order to test whether the results presented in this section could be artifacts of systematic effects that we failed to take into account, we sampled randomly from the posterior distributions of the parameters of the broken power-law model from all bursts where the null hypothesis was rejected. We then sampled the distributions of T90 and fluence from the observed ensemble of bursts and created fake light curves with the power spectral properties of the broken power-law model with the randomly sampled parameter sets and combinations of burst T90 and burst fluence taken from the real sample. We simulated light curves following Timmer & Koenig (1995) and included Poisson statistics in the simulated light curves. For 1000 such fake bursts, we fit the periodograms and performed model selection between the power-law and broken power-law models. A total of 707 periodograms simulated from the broken power-law model were actually adequately fit with a simple power law instead, whereas only 293 bursts rejected the null hypothesis. This indicates that the fit is strongly dependent on burst duration and fluence. Especially for short bursts, where there are few frequencies below 100 Hz, there may not be enough data points to require a more complex model. The power-law and broken power-law samples are well separated in terms of both burst T90 and burst fluence. A two-sided K-S test reveals a separation in burst duration, p = 1.0 × 10⁻⁴⁰, much stronger than observed in the data. The separation between the two samples in terms of the fluence distribution is comparable to that observed from the sample, p = 2.2 × 10⁻¹⁰. This indicates the possibility that separating bursts by their preferred broadband model may not be meaningful: perhaps all bursts follow a broadband noise process that is closer to a broken power law than a simple power law, but for short bursts, there are not enough data points at low frequencies to confidently reject a simple power-law model.

We also tested for the presence of a correlation between burst T90 and power-law index for the sample of fake periodograms that accepted a power-law model. While there still seems to be a correlation (Spearman rank coefficient p-value for no correlation: p = 1.35 × 10⁻¹⁰), this correlation is completely flat, with a power-law slope of α = 1.77 × 10⁻¹⁰, centered around a mean power-law index of 〈γ〉 = 2.56.

We have searched both individual bursts and averaged periodograms from samples of bursts for QPOs. Our analysis is the most precise to date for fast transients while taking into account the effects of aperiodic variability, but it is also conservative: at low frequencies, real quasi-periodic features may be missed because we assume that the burst is a purely stochastic process, when it is not. Additionally, we do not model several effects that can significantly affect the outcome of a QPO search. Dead time can significantly affect especially the bright bursts and thus render inferences invalid even at high frequencies. At low frequencies, we have demonstrated that some of the power concentrated in broad bumps, which can be classified as QPOs with broad widths by the detection algorithm, can easily be modeled with a simple estimate of a burst envelope consisting of two functions with a faster rise than decay. We must hence be very careful when interpreting signals at frequencies comparable to the lowest QPO frequencies seen in the giant flares (18 Hz and 36 Hz).

5.1. Individual Bursts

5.2. Averaged Bursts

5.3. Aperiodic Variability

The QPOs in the tail of the giant flare from SGR 1900+14 lie in the range 28–155 Hz (Strohmayer & Watts 2005). For the tail of the giant flare from SGR 1806−20, there are several QPOs in the range 18–150 Hz, and two isolated higher frequency signals at 625 Hz and 1840 Hz (Israel et al. 2005; Watts & Strohmayer 2006; Strohmayer & Watts 2006). Widths (FWHM) are in the range 1–20 Hz.

The most plausible explanation advanced for the giant flare QPOs is that they represent global seismic oscillations of the star, and it was swiftly realized that this would be a novel means of constraining not only the interior field strength (which is hard to measure directly) but also the dense matter equation of state (Samuelsson & Andersson 2007; Watts & Reddy 2007). The question of mode identification is therefore crucial.

In the original discovery papers, the QPOs were tentatively identified with torsional shear modes of the neutron star crust and torsional Alfvén modes of the highly magnetized fluid core. These identifications were based on the expected mode frequencies, which are set by both the size of the resonant volume and the relevant wave speed. For crustal shear modes, the appropriate speed is the shear speed v_s = (μ_s/ρ)^1/2, where μ_s is the shear modulus and ρ the density. The shear modulus is of the order of the Coulomb potential energy ∼Z²e²/r per unit volume r³, where r ∼ (ρ/Am_p)^−1/3 is the inter-ion spacing, while Z and A are the effective atomic number and mass number, respectively, of the ions in the crust. Using the shear modulus computed by Strohmayer et al. (1991) and scaling by typical values for the inner crust (Douchin & Haensel 2001), the shear velocity as shown by Piro (2005) is

$$\begin{eqnarray} v_s & = & 1.1 \times 10^8\ \mathrm{cm\ s^{-1}} \left(\frac{\rho }{10^{14} \ \mathrm{g \ cm}^{-3}}\right)^{1/6} \left(\frac{Z}{38}\right) \\ && \times \left(\frac{302}{A}\right)^{2/3} \left(\frac{1-X_n}{0.25}\right)^{2/3} \nonumber \,, \end{eqnarray} \tag{ 3 }$$

where X_n is the fraction of neutrons. This yields a rough estimate for the frequency for the fundamental crustal shear mode of ν ∼ v_s/2πR = 18 (10 km/R) Hz. Full-mode calculations find similar values, but with additional dependencies on the mass and radius of the star due to relativistic effects (see, e.g., Samuelsson & Andersson (2007)); it is this dependence that makes the modes potentially powerful diagnostics of the dense matter equation of state (Lattimer & Prakash 2007). Many of the lower QPO frequencies could be explained as angular harmonics with no radial nodes, while the two highest frequencies in the SGR 1806−20 giant flare were identified as radial overtones of these crustal modes.

For torsional Alfvén modes of the core, the appropriate wave speed is the Alfvén speed $v_A = B/\sqrt{4\pi \rho }$, where B is the magnetic field strength, giving

$$\begin{equation} v_A = 10^8 \ \mathrm{cm\ s^{-1}} \left(\frac{B}{10^{16}\ \mathrm{G}}\right) \left(\frac{10^{15} \ \mathrm{g\ cm}^{-3}}{\rho }\right)^{1/2}. \end{equation} \tag{ 4 }$$

This yields a very rough estimate for the frequency of the fundamental torsional Alfvén mode of ν ∼ v_A/4R = 25 (10 km/R) Hz (Thompson & Duncan 2001). Note, however, that the value of the field strength B in magnetar cores is highly uncertain, as is the appropriate value of the density ρ. In principle, only the charged component (∼5%–10% of the core mass) should participate in Alfvén oscillations, reducing ρ; however, there are mechanisms associated with superfluidity and superconductivity that can couple the charged and neutral components, leading to additional mass loading. As above, full-mode calculations that take into account relativistic effects lead to additional dependencies on neutron star mass and radius (see, e.g., Sotani et al. 2008). It should also be noted that the Alfvén modes constitute continua rather than a set of discrete frequencies, since the field lines within the core have a continuum of lengths. It has been suggested that the observed QPOs might be associated with a turning point in the continuum of Alfvén modes (Levin 2007; Sotani et al. 2008).

In fact, for a star with a magnetar strength field, crustal vibrations and core vibrations should couple together on very short timescales (Levin 2006, 2007). Considering them in isolation, as described above, is therefore not appropriate. The current viewpoint, based on more detailed modeling that takes into account the magnetic coupling between crust and core, is that the QPOs are in fact associated with global magneto-elastic axial (torsional) oscillations of the star (Glampedakis et al. 2006; Andersson et al. 2009; Steiner & Watts 2009; van Hoven & Levin 2011, 2012; Colaiuda & Kokkotas 2011, 2012; Gabler et al. 2012, 2013; Passamonti & Lander 2014, 2013; Glampedakis & Jones 2014).

Since magneto-elastic oscillations depend on the same physics described above, albeit now in a coupled system, they have frequencies in the same broad range as the simple estimates given above. The frequencies are set by many factors, including the dense matter equation of state (which sets mass and radius), field strength and geometry, superfluidity, superconductivity, and crust composition. Current magneto-elastic torsional oscillation models have had some success in explaining the presence of oscillations at frequencies of 155 Hz and below. However, they struggle to explain the presence of the highest frequency oscillations, which should damp very rapidly (van Hoven & Levin 2012; Gabler et al. 2012). Various solutions to this problem are under investigation, including coupling to polar modes (Lander et al. 2010; Lander & Jones 2011; Colaiuda & Kokkotas 2012), and resonances between crust and core that might develop as a result of superfluid effects (Gabler et al. 2013; Passamonti & Lander 2014). However, until these issues are resolved, precise identification of the giant flare frequencies with specific global magneto-elastic modes remains a challenge.

The detection of frequencies in the smaller flares provides an entirely new viewpoint on this very difficult theoretical problem. It is therefore important to compare the properties of the giant flare QPOs to those detected in this study. We begin with the frequencies detected by averaging together multiple bursts from highly active episodes, at 93 Hz and 127 Hz. These frequencies are in the range found in the giant flares (indeed the strongest frequency found in the SGR 1806−20 giant flare was at 92 Hz). The widths are also comparable to the range observed in the tails of the giant flares. It therefore seems plausible that they are instances of the same phenomenon. If these frequencies do indeed represent global magneto-elastic oscillations, the implication is that such vibrations are excited not only by giant flares but also by trains of shorter bursts. This is important information for future theoretical studies of mode excitation.

The 260 Hz signal identified in one of the individual bursts is rather different. It is found in a frequency range where no signals were found in the giant flares. It is much broader than any of the oscillations seen in the tails of the giant flares and has very high fractional amplitude. Whether it is the same phenomenon as was observed in the giant flares is therefore far from clear. If it is the same phenomenon, and we are seeing a magneto-elastic oscillation mode, then a detection in this frequency range would be valuable. Magneto-elastic oscillation models, as outlined above, have difficulty in explaining the lifetime of higher frequency signals in the giant flares. This detection, with a much lower coherence, in a burst whose duration is comparable to the predicted lifetimes, provides a fresh perspective on this problem.

We may, however, be seeing something quite different. The giant flares consist of a short impulsive spike, followed by a long decaying tail as a trapped pair-plasma fireball slowly evaporates. In the smaller bursts, it is not clear whether a fireball forms: what we see may be more analogous to the impulsive spike of the giant flares. The variability that we have found in the short bursts in SGR 1550 (particularly the 260 Hz signal that appears to differ in properties) may instead be associated directly with the burst trigger process, be that a magnetospheric instability or the yielding of the crust. It is worth noting that there are tentative claims of variability at 43 Hz in the first 200 ms of the 1979 giant flare from SGR 0526−66 (Barat et al. 1983) and at 50 Hz in the first 500 ms of the SGR 1806−20 giant flare (Terasawa et al. 2006; Geotail paper). However, timing analysis of the peaks of the giant flares is heavily affected by dead time and saturation. In this respect the smaller flares, which are typically not saturated, may be more useful despite the lower count rates. However, the variability that would be expected in the initial trigger and yielding phase of a magnetar burst has not been studied in detail. There are nonetheless plausible mechanisms that might lead to quasi-periodic behavior.

If the burst trigger is magnetospheric, there may arise via variability may rise via plasma instabilities associated with magnetic reconnection (see, e.g., Kliem et al. 2000). If instead the trigger is crustal yielding, local effects and resonances may be significant. It is not clear, for example, whether locally excited shear waves would immediately couple to the entire crust (and from there to the core) rather than resonating, with low Q-value, in a smaller cavity that is temporarily coupled very poorly to the rest of the crust during the yielding process. Such local resonances would have frequencies quite different from those of global magneto-elastic oscillations. More detailed theoretical studies of the trigger process will be required to resolve both this question and the length of time required to establish global modes of any kind. However, the possibility that the 260 Hz signal is a new and direct signature of the trigger process is an exciting one.

We have searched 263 bursts from SGR J1550−5418 for QPOs. We find one candidate QPO in the individual burst searches. The signal is broad, but highly significant, and not close to any frequency observed in the giant flares. It is unclear whether this signal could come from the same phenomenon as the QPOs observed in the giant flares, or whether it may be associated directly with the burst trigger process.

Searching averaged periodograms reveals a significant signal at ∼10 Hz in an averaged periodogram of all bursts with durations T90 > 0.5. This signal is comparable to the characteristic duration of a magnetar burst T90_max = 0.1 s, but may be due to an inadequate broadband model for the periodogram. We find evidence for QPOs in periodograms averaging bursts from individual triggers, which are unlikely due to effects of averaging together bursts with vastly different timing properties. Two of these signals are located at 93 Hz, where QPOs in the giant flares have been observed. The third is at a higher frequency, 127 Hz. We consider these signals to be the best candidates of neutron star oscillations from short magnetar bursts to date. All three signals can be interpreted in the framework of magneto-elastic oscillations invoked to explain QPOs in magnetar giant flares. The possibility that not only giant flares but also smaller flares may excite these oscillations also provides an important new piece of information for future theoretical studies of mode excitation. For averaged periodograms, we can put constraints on weak signals that could have been there and would likely have been missed by our methods. For all but the lowest frequencies, our sensitivity to QPOs is lower than the observed fractional rms amplitudes in the giant flares. This is especially prominent for the high-frequency QPOs observed at 625 Hz and higher. We thus conclude that except for the signals at 93 Hz and 127 Hz, there are no giant-flare-like QPOs in this sample of small bursts.

Here we also characterized overall burst variability for the first time. We find a correlation between power-law index and burst duration. This implies that longer SGR bursts are variable over a broader range of timescales than short bursts and are not simply longer versions of the short bursts. Further work is required to disentangle overall variability in magnetar bursts. This is unlikely to be possible with Fourier methods, but would be very rewarding in terms of both understanding emission mechanisms and untangling possible QPO signals from the overall burst morphology.

The authors thank Phil Uttley and Yuri Levin for useful discussions and the referee for helpful suggestions. D.H., C.D., and A.L.W. acknowledge support from a Netherlands Organization for Scientific Research (NWO) Vidi Fellowship (PI: A. Watts). C.K. was partially supported by NASA grant NNH07ZDA001-GLAST. This publication is part of the GBM/Magnetar Key Project (NASA grant NNH07ZDA001-GLAST; PI: C. Kouveliotou). A.J.v.d.H. acknowledges support from the European Research Council via Advanced Investigator Grant No. 247295 (PI: R.A.M.J. Wijers).