Stacking gravitational wave signals from soft gamma repeater bursts

Soft gamma repeaters (SGRs) have unique properties that make them intriguing targets for gravitational wave (GW) searches. They are nearby, their burst emission mechanism may involve neutron star crust fractures and excitation of quasinormal modes, and they burst repeatedly and sometimes spectacularly. A recent LIGO search for transient GW from these sources placed upper limits on a set of almost 200 individual SGR bursts. These limits were within the theoretically predicted range of some models. We present a new search strategy which builds upon the method used there by “stacking” potential GW signals from multiple SGR bursts. We assume that variation in the time difference between burst electromagnetic emission and burst GW emission is small relative to the GW signal duration, and we time-align GW excess power time-frequency tilings containing individual burst triggers to their corresponding electromagnetic emissions. Using Monte Carlo simulations, we confirm that gains in GW energy sensitivity of $N^{1/2}$ are possible, where $N$ is the number of stacked SGR bursts. Estimated sensitivities for a mock search for gravitational waves from the 2006 March 29 storm from SGR 1900+14 are also presented, for two GW emission models, “fluence-weighted” and “flat” (unweighted).

Figure 1

BAT light curve of the SGR 1900+14 storm event, 1 ms bins, 15–100 keV. The light curve shows the BAT event data, from sequence 00203127000, from approximately 20 s after the start of the sequence to its end. The data are publicly available [42]. Several of the largest bursts, durations greater than 500 ms, are considered intermediate flares. The $x$-axis times are relative to 2006 March 29 $02:53:09.9\pm 0.5\mathrm{s}$ UT at the satellite. Times used in the analysis itself have been corrected for time-of-flight delays to the geocenter, which change by 0.35 ms over the $\sim 30.5\mathrm{s}$ course of the light curve.

Figure 2

Individual EM bursts inform the stacking of GW data. This figure suggests the stacking procedure and explicitly shows search time scales. The top four plots are EM light curve time series around individual bursts beginning at 2.0 s, 16.6 s, 18.3 s, and 22.3 s in Fig. 1. Simulated GW ringdowns in the fluence-weighted model are superposed. The bottom plot shows the EM time series simultaneously, and the sum of the hypothetical GW signals. The on-source region of $\pm 2\mathrm{s}$ is shaded. In a search, GW data corresponding to the EM time series are transformed to time-frequency power tilings before being added together and therefore there is no dependence on phase-coherence of GW signals in the analysis; this transformation is not illustrated.

Figure 3

Diagram of the $T$-stack version of the stack-a-flare pipeline. The $T$-stack pipeline has a thin layer added before the flare pipeline in which gravitational wave data time series containing SGR burst event triggers are aligned on the trigger times and added together. These stacked time series are made for each detector and then run through the flare pipeline as normal.

Figure 4

Diagram of the $P$-stack version of the stack-a-flare pipeline. The $P$-stack pipeline has a thin layer added after the flare pipeline in which gravitational wave data significance tilings containing SGR burst event triggers are aligned on the trigger times, weighted according to astrophysical factor, and added together. Stacked significance tilings can then be run through the flare pipeline clustering algorithm.

Figure 5

Example efficiency curve generated for the Monte Carlo experiment investigating stack-a-flare sensitivity vs $N$. This example curve was made with the $T$-stack pipeline using compound injections of $N=200$ 1590 Hz circularly polarized ringdowns injected into Gaussian noise with $\sigma =1$. Each efficiency curve was constructed using 20 amplitude scaling factors and 20 trials at each $h_{\mathrm{rss}}$ amplitude. Binomial error bars are given by $\sqrt{r(1-r)/M}$ where $r$ is the fraction detected at a given injection amplitude and $M$ is the number of injections in the Monte Carlo.

Figure 6

$T$-stack and $P$-stack sensitivity dependence on $N$ at 90% detection efficiency, for 1590 Hz $\tau =200\mathrm{ms}$ ringdowns in Gaussian noise with $\sigma =1$. The $T$-stack pipeline shows a sensitivity dependence of nearly $N^{1/2}$ and the results for the $P$-stack pipeline show a sensitivity dependence of nearly $N^{1/4}$. All fits excluded the point $N=1$. Results at 50% detection efficiency look similar.

Figure 7

$T$-stack and $P$-stack sensitivity dependence on $N$, at 90% detection efficiency, for 100–1000 Hz 100 ms duration white noise bursts in Gaussian noise with $\sigma =1$. The results for the $T$-stack pipeline show a sensitivity dependence of nearly $N^0$ (flat dependence), and the results for the $P$-stack pipeline show a sensitivity dependence of nearly $N^{1/4}$, as in the coherent ringdown case. All fits exclude the point $N=1$. Results at 50% detection efficiency look similar.

Figure 8

$T$-stack and $P$-stack sensitivity versus timing error, for 1090 Hz $\tau =200\mathrm{ms}$ circularly polarized RD, $N=20$, for $h_{\mathrm{rss}}^{90\%}$. $T$-stack is more sensitive for small timing errors but degrades. The crossover point is noted; for $T$-stack to be effective at 1090 Hz, apparently timing error must be $\lesssim 100\mu \mathrm{s}$ at one-sigma. $T$-stack results level off at high timing errors (greater than $\sim 2\times {10}^{-4}$, or $\sim 90\mathrm{degrees}$ of phase) because the Monte Carlo effectively randomizes the phases of the stacked signals. The shaded region shows the approximate range of timing uncertainties, shown in Table . Results at 50% detection efficiency look similar.

Figure 9

$T$-stack and $P$-stack sensitivity versus timing error, for 2590 Hz $\tau =200\mathrm{ms}$ circularly polarized RD, $N=20$, for $h_{\mathrm{rss}}^{90\%}$. $T$-stack is more sensitive for small timing errors, but degrades. The crossover point is noted; for $T$-stack to be effective at 2590 Hz, apparently, timing error must be $\lesssim 50\mu \mathrm{s}$ at one-sigma. $T$-stack results level off at high timing errors (greater than $\sim 1\times {10}^{-4}$, or $\sim 90\mathrm{degrees}$ of phase) because the Monte Carlo effectively randomizes the phases of the stacked signals. The shaded region shows the approximate range of timing uncertainties, shown in Table . Results at 50% detection efficiency look similar.

Figure 10

SGR 1900+14 storm light curve with 1 ms bins in the (15–100) keV band. Bottom plot shows a detail. Burst start times are estimated by fitting the steeply rising burst edges; EM fluences are estimated by integrating light curve area under each burst. A 30-bin running average is shown in addition to the raw light curve. Solid lines are linear fits to rising edges; the boundaries of rising edges were found by examining the first derivatives in the neighborhoods of the peak locations. Crosses mark burst peaks and intersections of the rising edge fits extrapolated to a linear fit of the noise floor measured in a quiescent period of data in the 50 s BAT sequence before the start of the storm. The one-sigma timing uncertainty averaged over all measurements is 2.9 ms. $X$-axis times are relative to 2006 March 29 02:53:09.9 UT at the Swift satellite.

Figure 11

Histogram of one-sigma timing uncertainties associated with the rise start time of each burst. These uncertainties were estimated from linear fits of the rising edges of bursts in the light curve assuming errors in the light curve data are independent normal with constant variance. The smallest uncertainty is $1.6\times {10}^{-04}\mathrm{s}$, the largest is 18 ms, and the mean is 2.8 ms. The largest uncertainty is due to a burst at $\sim 4.4\mathrm{s}$ which rises out of another large burst, far above the noise floor. The other large uncertainties are due to bursts with rippling rising edges which may be due to additional injections of energy.

Figure 12

SGR 1900+14 storm light curve, focusing on fluence (integrated counts) estimates and time scales. The $N=11$ flat model is illustrated here. In addition to features common to Fig. 10, alternating red and cyan highlights indicate integration bounds for fluence estimation of the 11 most energetic bursts. Fluence estimation is one of the goals of the light curve processing algorithm. We show $\tau =200\mathrm{ms}$ ringdown envelopes, aligned with the nominal beginning of the 11 brightest electromagnetic bursts. A small constant time offset between this nominal time and gravitational wave emission would be immaterial to the search, so long as the time offset were small relative to the on-source region half-duration of 2 s. The $\pm 2\mathrm{s}$ on-source regions, which often overlap even in this $N=11$ model, are illustrated by the vertical lines and shading.

Figure 13

Histogram of integrated counts in bursts in the SGR 1900+14 storm. This quantity is roughly proportional to fluence. The histogram shows a separation between the 11 most electromagnetically energetic bursts and the rest. This separation determined the break point for the $N=11$ flat GW emission model.

Figure 14

Top plot: cumulative integrated counts. The vertical line marks the contribution from the 11th most energetic peak. 89% of the total integrated counts in the 42 storm bursts measured is contained in the most energetic 11 bursts. The 12th most energetic burst contributes an additional 1% to the total estimated from the 42 bursts. Bottom plot: fractional change in the cumulative integrated counts. The 16th burst contributes less than 1% to the running total. The 37th burst contributes less than 0.1% to the running total.

Figure 15

Example cumulative histograms showing false alarm rates versus analysis event loudness for the background (blue) and the stacked on-source region (five red points). There are three such plots for each emission model, one per search band. Since the stacked live time is 4 s here, the loudest on-source event occurs once per 4 s, and is plotted at a $y$ value of 0.25 Hz. We can estimate the FAR of this loudest on-source analysis from the background. Error bars are Poissonian at 90% confidence.

Figure 16

Stack-a-flare simulated data energy sensitivity estimates at 10 kpc, for the 29 March 2006 storm from SGR 1900+14, for the flat and fluence-weighted models. Uncertainty estimates have been folded in, as tabulated in Table . Vertical lines indicate boundaries of the three distinct search frequency bands. Crosses and circles indicate linearly and circularly polarized RDs, respectively. Triangles and squares represent 11 ms and 100 ms band- and time-limited WNBs, respectively. Symbols are placed at the waveform central frequency. These sensitivity estimates reflect the noise curves of the detectors.