Broad band turbulent spectra in gamma-ray burst light curves
Broad band power density spectra offer a window to understanding turbulent behavior in the emission mechanism and, at the highest frequencies, in the putative inner engines powering long GRBs. We describe a chirp search method which steps aside Fourier analysis for signal detection in the Poisson noise-dominated 2 kHz sampled BeppoSAX light curves. An efficient numerical implementation is described in operations, where is the number of chirp templates and is the length of the light curve time series, suited for embarrassingly parallel processing. For detection of individual chirps of duration s, the method is one order of magnitude more sensitive in SNR than Fourier analysis. The Fourier-chirp spectra of GRB 010408 and GRB 970816 show a continuation of the spectral slope up to 1 kHz of turbulence identified in low frequency Fourier analysis. The same continuation is observed in an ensemble averaged spectrum of 40 bright long GRBs. An outlook on a similar analysis of upcoming gravitational wave data is included.
I 1. Introduction
High frequency power density spectra up to 1 kHz offers a window to the inner engines of long GRBs, that may harbor a rapidly rotating (proto-)neutron star (PNS) or a black hole-accretion disk or a torus system (BHS; e.g. pir99 (); pir04 ()). Such high frequency window may be probed by analyzing the 2 kHz light curves in the BeppoSAX catalogue fro09 () or the upcoming strain amplitude data from gravitational wave detectors LIGO-Virgo, KAGRA and the Einstein Telescope bar99 (); arc04 (); kur10 (); hil08 ().
Given the limited collector area of today’s gamma-ray satellites, high frequency light curves of GRBs are typically Poisson noise dominated, flattening their Fourier spectra above at most tens of Hz. Low-frequency Fourier analysis reveals a Kolmogorov spectrum in GRB light curves bel98 (); bel00 (); gui12 (); dic13 (), that is expected to continue smoothly to high frequencies. A broad-band turbulent spectrum from the gamma-ray emission process may hereby present a new baseline in searches for high frequency modulations by the central PNS or BHS.
Here, we describe a method to extend the turbulent spectra to high frequencies in the Poisson noise of BeppoSAX light curves of long GRBs. Turbulence produces phase-coherent intermittencies on short to intermediate time scales, that covers an extended bandwidth in frequency space. To search for turbulence in a Poisson noise dominated signal, we therefore set out to step aside Fourier analysis, since it focuses on phase coherence across a narrow bandwidth. We shall apply matched filtering with chirps whose frequencies increase or decrease exponentially in time, representing phase coherence across a finite bandwidth.
The computational effort is operations, consistent with the Fast Fourier Transform, where denotes the number of samples in the data time series and the number of chirp templates. This limit appears to be competitive to other approaches oto10 (); fis11 ().
To illustrate our method, we report on spectra of bright 2 kHz sampled light curves of long GRBs in the BeppoSAX catalogue fro09 ().
The chirp templates are described in §2 and an efficient numerical implementation of the matched filtering algorithm is given in §3. Chirp detection sensitivity is analyzed in §4. We apply the method to extract broad band spectra from a sample of BeppoSAX light curves in §5,6, to derive a continuation of Fourier spectra. An outlook on further applications is briefly described in §7.
Ii 2. Chirp templates
Chirps are transients described by a base frequency and a frequency rate of change. They are different from quasi-periodic oscillations (QPOs), that pertain to frequencies meandering about a steady mean. QPOs are notoriously absent in GRB light curves cen10 (), whence they will not be considered here.
Given the limitations of Fourier analysis to extract broad band turbulent spectra in noisy data, we here consider a search by matched filtering for time coherent features across a finite frequency bandwidth. Of particular interest are chirps, here with an exponential change in frequency as a function of time. They can be extracted by time slicing a single long duration chirp van11 () into subintervals of duration much shorter than the duration of a GRB. Starting from frequency and decaying to a late time, asymptotic frequency , the frequency evolution is essentially exponential in time,
where is a dimensionless scale to parametrize the time scale of change in chirp frequency.
The choice of is used to search for phase coherence over a time scale , such that , where denotes the sampling time interval or the bin size of integration of the data. Accordingly, we consider time intervals
in our time-slicing procedure. Fig. 1 illustrates slicing of a model chirp with s into one second chirp templates
In steady state, the statistical properties of, e.g., a velocity field in turbulent motion are the same when viewed forwards and backwards in time, whereby the probabilities for detecting positive and negative chirps are similar for intermediate durations . We shall therefore employ difference chirps, obtained as the difference of chirps considered forwards and backwards in time. If is a chirp template extracted from (5), we consider
In using (8), we save a factor of two in computational cost, allowing for detections of positive or negative chirps in one calculation. Because the cross-correlation between chirps forward and backwards in time is small, (8) can be used efficiently with negligible loss in sensitivity over performing matched filtering twice, using positive and negative chirps separately in each run.
Iii 3. Efficient matched filtering
To develop a search with maximal sensitivity to features with frequencies changing in time, we set out to employ matched filtering. Matched filtering obtains the highest possible sensitivity for phase-coherent features, upon including sufficiently many templates to cover the full range of possible signals and their phase-coherent behavior. For this reason, efficient numerical implementation of method matched filtering is important.
For a time series and chirp template , let and by subtracting the respective mean values and , and consider their cross correlation
Using the Fourier transform of a function ,
is efficiently calculated as the inverse Fourier transform of the product of the complex conjugate of the transform of and of .
In matched filtering, the potential significance of a chirp is identified by the normalized cross correlation of and (i.e., the Pearson coefficient)
where is the norm of .
For discrete series of samples at , consider and , where denotes the number of slices illustrated in Fig. 1. Since (5) is bilinear, the sample correlation coefficient (SCC) obtained from on and is readily calculated by the Fast Fourier Transform (FFT) in operations. However, the normalized sample correlation coefficient obtained from and is nonlinear in and due to normalizations . Direct evaluation of these nonlinear expressions is prohibitively expensive when the number of chirps is large.
To be precise, consider the Pearson coefficient between given by a section ), of the time series and a template ,
Exact normalization in (9) ensures it to be the cosine between and , satisfying
The unnormalized cross correlation
can be evaluated in operations using the Fast Fourier Transform (FFT), which is is cost effective compared to (9) when . We now consider block wise over data slices of duration , i.e., the discretized intervals
Normalizing as an array over each slice requires normalizations,
The define a staircase as a function of , where , giving a block wise normalized SSC
Block wise normalization is particularly opportune when shows variations in dispersion no faster than the intermediate time scale . If so, our block normalized SCC applies at the operational cost of FFT, which is below the cost of true SCC’s in (9) whenever .
When is relatively large, e.g., with , typically shows a near-Gaussian distribution by the central limit theorem. In what follows, we shall use a further normalization by the variance of the cross correlations (11) in each time slice, i.e.,
hereby differs from only by a constant factor , whose Probability Density Function (PDF) approaches a truncated Gaussian of unit variance. The truncation is a function of both the total number of trials and the potential for a chirp being present in the data.
Fig. 2 illustrates the numerical implementation.
Iv 4. Chirp detection sensitivity
We consider quantifying the sensitivity of matched filtering relative to that obtained by Fourier analysis. To this end, we perform the injection experiment
on the BeppoSAX light curve of GRB 010408 by the light curve of a chirp for various amplitudes .
A chirp search by matched filtering obtains the time series following (15). We consider detection by
in light of the approximately Gaussian PDF of the , truncated by the finite number of trials in each template search.
To compare (17) with Fourier analysis, we calculate the spectrum by the Welch method wel67 (); coo70 (); pre02 () using a partition in sub-windows of length with a Hz frequency resolution. The spectrum is calculated as an average over periodograms from intervals of length with 50% overlap. Each periodogram is obtained with a Welch window, i.e., as the Fourier transform of by FFT, where , where with . By construction of the Welch method, fluctuations in the resulting power spectrum are approximately Gaussian, as a distribution from averaging periodograms.
For each , detection is expressed, similarly to (17), by peak values
relative to the asymptotically flat Poisson dominated spectrum in terms of
where and denote the mean and standard deviation of the in the Poisson dominated tail of the Fourier spectrum.
As a control, consider a chirp search in random data. In this event, and have expectation values and , respectively, that derive from the expectation value of the truncation in the distribution of trials of a variable taken from a Gaussian distribution with unit variance. Here, satisfies
with for and , for , where denotes the complementary error function.
Fig. 3 shows and distributions computed numerically from maxima in trial samples of size from a Gaussian distribution with unit variance. Their mean is determined by according to (20) and the resulting distributions have positive skewness. The distribution shown are themselves truncated according to (20) upon substituting for .
Figs. 4 and 5 show the results of our injection experiment for various values of obtained in Fourier analysis and, respectively, matched filtering.
Applied to , (20) implies a base level when in considering 8 seconds of the 2 kHz light curves in the BeppoSAX catalogue. Under the null hypothesis of no signal present, has a probability of occurrence . An excess is a false positive with or denotes the presence of a signal. For ( s), a detection corresponds to , indicated by the dot-dashed line in Fig. 5. obtains similarly from giving a threshold .
For a detection, the critical value for matched filtering is a factor of 5.3 smaller than for Fourier analysis using the Welch method. An additional moving average of the Fourier spectrum shows some improvement in its sensitivity, leaving a gain by matched filtering by about one order of magnitude in SNR.
The sensitivity of matched filtering shown in Fig. 5 is consistent with the theoretical sensitivity limit of , where in the case at hand for s. The observed value for a 1- excess in above background is indeed close to the anticipated value .
V 5. Chirp spectra
We observe that both scales linearly with across an appreciable range beyond , below which it assumes the constant background value set by the number of trials . Upon subtracting , we thus obtain a detection method with a linear response to small amplitude signals. In general, is a function of frequency, which poses the question on devising a suitable control.
We express the spectra in terms of a strain given by the square root of the PDS,
Here, denotes the results of chirp analysis of control light curves and is the bandwidth of a chirp templates about the frequency , where is a coefficient that depends on the choice of chirp duration. The are calculated from the Fourier coefficients , that have a noise dominated high frequency tail with standard deviation and mean .
A chirp search is not a linear transform in the sense of Fourier analysis. A chirp search seeks a best-fit chirp to the data, by varying frequency and frequency rate of change. Different chirp templates are hereby linearly dependent at high resolutions even though they may retain finite cross-correlations, as opposed to working with basis functions that satisfy exact linear independence. However, this distinction is immaterial in calculating spectra.
For Gaussian additive noise, such as in the high frequency, shot-noise dominated region of strain amplitude noise in gravitational wave detectors, light curves obtained by time randomization or produced by a random number generator will be effective as a control . By whitening in (15), these two alternatives give essentially the same results.
For the Poisson noise in the 2 kHz BeppoSAX light curves, whose average photon counts are of the order of unity per 0.5 ms bin, we propose as a control a synthetic Poisson noise light curve, that shares the same smoothed light curve as the original. A control of this type accurately captures the secular variation of the variance in the noise with (smoothed) amplitude, illustrated in Fig. 6.
We are now in a position to apply our method to two bright long GRBs from the BeppoSAX catalogue. Matched filtering calculations are performed with a log-unfiorm distribution in frequency and frequency rate of change over a total of 5.76 million templates. For control, we use synthetic Poisson noise light curves about smoothed light curves following a low-pass filter at 2 Hz.
Fig. 7 shows a blended Fourier-chirp spectrum up to about the maximal frequency of 1000 Hz set by the sample rate of 2 kHz. Here, the low frequency spectrum is computed by Fourier analysis and the high frequency spectrum by our chirp search method. Included is a linear extrapolation of the low frequency spectrum, to highlight a common spectral slope in the low and high frequency spectra, here obtained independently by two completely different methods.
Vi 6. Ensemble spectrum of bright GRBs
We next consider a sample of 72 bright events in the BeppoSAX catalogue (Fig. 8). We select a subsample of 40 events with a pronounced autocorrelation in their 2 kHz light curves (“red,” Fig. 9) for extracting an ensemble averaged Fourier-chirp spectrum.
The Swift catalogue of long GRBs shows no correlation between the observed durations redshift, shown in Fig. 10. The spread in the observed durations , therefore, is, essentially intrinsic to the source.
Based on van11b (), we consider long GRBs to be produced by black hole-disk or torus systems (BHS), rather than a (proto-)neutron star, based on two hyper-energetic GRB-supernovae with an output exceeding the maximal spin energy of the latter. In a BHS, can be identified with the lifetime of black hole spin, whereby for a black hole mass van01 (). The spectrum of turbulence and intermittencies in the surrounding accretion disk or torus scales likewise with . If correlated to the wind from the disk or torus, the spectrum of the turbulent outflow is normalizable by multiplication by , i.e., by .
Fig. 11 shows the ensemble average of the normalized spectra, plotted as a function of normalized frequency in the co-moving frame of reference, assuming a fiducial redshift similar to the mean of in the Swift sample shown in Fig. 10.
A detailed consideration of alternative chirps shows the following. Chirps with constant amplitude produce slightly more scatter than those obtained from time slicing shown in Fig. 1. Chirps with s ( s) show considerably more (less) scatter in the ensemble average. Our choice of appears to provide a compromise between scatter and frequency coverage in extending the slope of the turbulent spectrum.
Vii 7. Conclusions
Turbulent spectra in low frequency in Fourier analysis of long GRBs are found to have a continuation to high frequencies, here found in two relatively bright GRBs obtained in a broad band chirp search by matched filtering. Matched filtering theoretically obtains maximal sensitivity for a detection, provided that the template bank is sufficiently dense and broad to capture the signal of interest. To capture turbulence, we here employ chirps with frequencies varying slowly in time following exponential decay or growth. Fig. 7 shows that Poisson noise is hereby effectively circumvented, upon using a control that shares the secular evolution of variance with amplitude of Poisson noise in the 2 kHz BeppoSAX light curves (6).
Fig. 7 shows that extraction of high frequency spectra are quite noisy due to the strong Poisson noise in the BeppoSAX light curves in light of the small number photon counts in each 0.5 ms bin. On this basis and our limited chirp parameter scan, e.g., using only, there is no conclusive evidence for the presence or absence of pronounced transient chirps distinct from those arising from turbulence. Extensive searches for transient chirps of different durations await a future investigation.
Our extension of the turbulent spectrum in the blended Fourier-chirp spectrum can serve as a new base line in searches for high frequency transient features. To this end, we consider the smoothed ensemble averaged spectrum of Fig. 8. It may serve as a reference in searches for bumps at high frequency, e.g., around the de-redshifted frequency of 1 kHz. Detection of a bump would reveal the presence a PNS with misaligned axis of angular momentum and magnetic field, representative for the birth of a new pulsar. The same would be absent in case of rapidly rotating black hole, whose magnetic moment and angular momentum are perfectly aligned by Carter’s theorem car68 (). Based on the present chirp search over 5.76 million templates, no such bump is found.
Figs. 7 and 11 demonstrate high frequency analysis as a new probe of the physics of the gamma-ray emission mechanism, that includes a potentially powerful window to intermittencies in the GRB inner engine even in light of exceedingly small photon counts. The ensemble of 40 bursts used in Fig. 11 represents bright events with a pronounced autocorrelation (“red” events with mean photon counts of 1.2569 per bin in the brightest channel) selected out of an initial sample of 72 bright events in the BeppoSAX catalogue, the remaining 32 (“white” events with mean photon counts of 0.5936 per bin in the brightest channel) lacking any perceptible autocorrelation with photon counts lower by a factor of about two. Thus, future gamma-ray missions with larger photon yields promise to greatly facilitate high frequency analysis, by improving signal-to-noise ratios and enlarging the sample of red events.
Chirp searches can also be applied the strain amplitude data from upcoming advanced gravitational wave detectors LIGO-Virgo and KAGRA, by changing control to, e.g., time randomized data. The proposed Fourier-chirp spectra can be extracted to search for gravitational wave signatures of possibly forced turbulence van01 () in high density matter in the inner disk or torus around black holes, long lasting over up to tens of seconds and possibly accompanied by pronounced transient chirps van11 (). Given the limited sensitivity range of these detectors, of interest are LGRBs and hyper-energetic core-collapse supernovae in the Local Universe. Core-collapse supernovae may be found in nearby galaxies such as M51 (hosting SN1994i, SN2005cs, SN 2011dh) and possibly M82 mux10 () with event rates over one per decade in each.
Acknowledgment. The BeppoSAX mission was an effort of the Italian Space Agency ASI with participation of The Netherlands Space Agency NIVR. Computations in this research were supported in part by the National Science Foundation through TeraGrid (now XSEDE) resources provided by Purdue University under grant number TG-DMS100033. Some calculations were performed at CAC/KIAS and KISTI. F. F. and C. G. acknowledge financial support from Italian Ministry of Education, University and Research through the PRIN-MIUR 2009 project on Gamma Ray Bursts (Prot. 2009 ERC3HT).
- (1) B. Barish and R. Weiss, Phys. Today 52, 44 (1999)
- (2) F. Arcese et al., Classical Quantum Gravity 21, S385 (2004)
- (3) K. Kuroda et al. (LCGT Collaboration), Class. Quantum Gravity 27, 084004 (2010)
- (4) C. Guidorzi, M. Margutti, L. Amati, et al., MNRAS, 422, 1785 (2012)
- (5) S. Hild, S. Chelkowski, and A. Friese, arXiv:0810.0604
- (6) M.H.P.M. van Putten, N. Kanda, H. Tagoshi, D. Tatsumi, F. Masa-Katsu, & M. Della Valle, Phys. Rev. D, 83, 044046 (2011)
- (7) T. Piran, PhR, 314, 575 (1999)
- (8) T. Piran, RvMP, 76, 1143 (2004)
- (9) A.M. Beloborodov, B.E. Stern, R., Svensson, ApJ, 508, L25 (1998)
- (10) A.M., Beloborodov, B.E. Stern, R. Svensson, ApJ, 535, 158 (2000)
- (11) S. Dichiara, C. Guidorzi, F. Frontera & L.A. Amati, ApJ, 777, 132
- (12) O. Korobkin, E.B. Abdikamalov, E. Schnetter, N. Stergioulas & B. Zink, PRD, 83, 043007 (2011)
- (13) J.M. OÕToole, M. Mesbah, & B. Boashash, IET Signal Process., 4, 428 (2010)
- (14) A. Fish, S. Guevich, R. Hadani, A. Sayeed & O. Schwarz, arXiv:1112.4883v1 (2011)
- (15) Cenko, S.B., et al., 2010, ApJ, 140, 224
- (16) F. Frontera, C. Guidorzi, E. Montanari, et al., ApJ Suppl., 180, 192 (2009)
- (17) P.D. Welch, IEEE Trans. Audio Electroacoustics, AU-15, 70 (1967)
- (18) W. H. Press, S. A. Teukolsky, W. T. Vetterling & B. P. Flannery, Numerical recipes in C, The art of scientific computing (Cambridge: Cambridge University Press, 2002), §13.7
- (19) Cooley, J.W., Lewis, P.A.W., & Welch, P.D., J. Sound Vob., 12, 339 (1970)
- (20) NASA, HEASARC, http://swift.gsfc.nasa.gov/archive/grbtable
- (21) van Putten, M.H.P.M., Della Valle, M., & Levinson, A., 2011, A&A, 535, L6
- (22) van Putten, M.H.P.M., & Ostriker, E.C., ApJ, 552, L31
- (23) B. Carter, Phys. Rev., 174 1559 (1968)
- (24) M.H.P.M. van Putten, 284, 115 (1999); ibid. Phys. Rev. Lett., 84, 091101 (2001); ibid. ApJ, 575, L71 (2002)
- (25) T.W.B. Muxlow, R.J. Beswick, S.T. Garringon, et al., MNRAS, 404, L109 (2010)