Distinguishing Time Clustering of Astrophysical Bursts

The time evolution of energetic astrophysical burst events carries critical information and clues to their nature. This can reveal exciting, extreme astrophysics such as complete stellar disruption, ultrahigh magnetic fields, accompanying neutrino bursts, etc. and high luminosities visible to great distances. Notable high luminosity examples include gamma ray bursts (GRB), fast radio bursts (FRB), and soft gamma repeaters (SGR). Important puzzles remain concerning their nature, involving aspects of stellar structure, accretion, jet production, and the circumobject medium.

Of particular interest are those events that recur, indicating that the process is not wholly disruptive, and possibly involves rotation or orbits; they also offer the possibility of observing the explosive event multiple times. Furthermore, if the repetition can be somewhat predicted, so that observations can be scheduled, this enables enhanced opportunities for understanding the burst mechanism and astrophysics. Orbits and rotation naturally impose periodic modulation on a variety of astrophysical signals, from occultations of stars by planets to pulsars to accretion disk phenomena. As a result there is a highly developed set of mathematical, statistical, and computer code tools to estimate the statistical likelihood of periodicity in a noise-limited or non-ideal sampling of data, and then determining the period and its uncertainty. This field has shown diverse evolution and activity from early Blackman-Tukey analysis [1] to “pulsar folding” (e.g. [2]) to Lomb-Scargle periodograms (see, e.g., [3] for a review of the Lomb-Scargle periodogram and comparison to other methods).

An example of greater complexity of phenomena and data came with the discovery and study of quasi-periodic oscillations – such an X-ray emission phenomenon linked to accretion appears not as a simple periodic signal but as a broad bump in frequency space. Our focus here is periodic windowed behavior (PWB), inspired by the recent discovery of such behavior in repeating FRB sources [5, 4]. For PWB, activity occurs only during periodically occurring windows – there is no activity in the gaps between the active windows, however not all active windows may show activity. The activity within a window may be random; there is no requirement or expectation that it will converge to a uniform profile (like a pulsar profile). Both the period and the active fraction (e.g. related to duty cycle of the energetic astrophysical process) are of interest.

Radio telescopes observing repeating FRB sources report millisecond duration bursts from the same source that may have time spacings from milliseconds to days during continuous observations, but intensive monitoring campaigns may observe no bursts for $\sim 180$ days [4, 5]. Initially, it was proposed that the behavior of the best-known repeating FRB, the source of FRB 121102, may be modeled as a time-clustered Weibull distribution [6]. However, with a few years of data in hand, an unexpected, and stunning, result appeared: the bursts were not periodic, or simply clustered, but were observed only in“periodic activity windows” [4, 5]. Similar behavior was reported for the source of FRB 180916 [7].

More recently, PWB was reported in SGR1935+2154, but in soft gamma-ray bursts [8]. This object is particularly noteworthy as it became the only known source of FRBs within our Galaxy when two were detected during the same day in 2020 [9, 10]. This object is so much closer than any other known FRB source that one could hope it could be some kind of “Rosetta stone” for understanding FRBs. (Another claim of PWB in soft gamma bursts was given for SGR1806-20 [11], though it is less certain.)

Robust identification, and characterization, of PWB could shed light on astrophysical burst mechanisms and energetics. For example, a considerable number of theories have been proposed for the origin of FRBs [12]; documentation and analysis of PWB in FRBs would be very important in constraining these models and could be key to understanding the bursting nature. As we enter an age of big data time domain surveys, PWB could also be discovered in other astrophysical contexts, showing further value for improved ways of identifying and measuring this phenomenon.

In Section II we describe a method for identifying PWB, and in particular distinguishing it from more irregular time clustering of events. For a burst series suspected of having PWB, in Section III we present methods for determining the period and active fraction. We apply this to the data from SGR1935+2154 in Section IV, and discuss extensions and conclude in Section V.

Strictly periodic behavior can be identified through a large number of different methods; for astrophysical time series data one of the most widely used is the Lomb-Scargle periodogram [13, 14, 3], which works in frequency space, where the folding of the time series increases the signal. In periodic windowed behavior, activity occurs at some times (not necessarily deterministic) within windows, where it is the windows that recur at a regular period. Such behavior may be the result of, e.g., a periodic “shutter” modulating non-periodic emission, or the physical conditions necessary for an outburst may occur only periodically (but not guarantee a specific time for a burst), or see [15] and references therein for an FRB precessing beam model. In this case long-term folding of data will not necessarily converge to an average profile.

The active windows may be a significant fraction of the full period, and then often the frequency analysis methods are diluted, with aliasing of the signal, i.e. signal peaks are broadened. Moreover, they tend to ignore the “curious incident of the dog in the night time” [16]: that the dog didn’t bark (the gaps with no bursts) carries important information. Therefore we discuss a time domain method for assessing the burst data.

Working in the time domain, we initially examine the cumulative distribution function (CDF), since it does contain information on both events and gaps, and has useful statistical properties. From an ordered series of events at times $\{t_{i}\}$, we look at the distribution of burst spacings $t_{i}-t_{1}$ vs observation time. Normalizing both axes to run from 0 to 1 (at the end we can restore the scaling by the total observation time, to get answers in days), we have the cumulative fraction of events, between 0 and 100%, vs the time fraction, normalized by the length of the observation campaign¹¹1Some technical details: we take time 0 to be when the first burst is observed; any gap before this is uninformative since we have no way of knowing whether the object had ever burst before – basically we are interested in the repeats. The end time can be either the end of the observing or when the last burst is detected; this does not affect the results. We then normalize all burst time differences $t_{i}-t_{1}$ by $t_{\rm end}-t_{1}$ so that the $x$-axis of time fraction runs from [0,1].. Thus we have the CDF of a variable distributed on [0,1].

Figure 1 illustrates several CDF for various realized distributions with a mean of 128 events each. We compare the uniform random distribution to two clustered distributions – all without actual periodicity – to a PWB case. The distributions have both similarities and differences, and the eye can be fooled by random excursions into thinking it detects patterns, even periodic ones.

The uniform distribution is realized by generating event times in a uniform random fashion over the [0,1] interval. The Weibull distribution is often used in astrophysics (e.g.  [6, 5] for FRBs), and more widely to give events that are clustered together at either early times or late times (as in failure rapidly or due to aging). We show an example with scale parameter $\beta=0.4$ and shape parameter $k=1.5$, so the CDF is $1-\exp\{-(t/\beta)^{k}\}$. Events beyond the end of the unit interval are not yet observed and not selected. The clustering here is due solely to the distribution. A true clustering is implemented in the clustered distribution, where events are generated with a two point correlation function $\xi(r\equiv|t_{j}-t_{i}|)\sim r^{-\gamma}$. We construct this using the Soneira-Peebles approach [17, 18], using a multiplicity parameter $\eta=4$, scale parameter $\lambda=8$ (thus $\gamma=0.33$), and four levels, giving 256 points, from which we randomly select the desired number.

For the PWB distribution we take four periods of length 0.25 with active fraction 60%, i.e. an activity length of 0.15 and a gap length of 0.1. Within the active window the events are uniformly randomly distributed. While by eye one can discern gaps of no activity in the PWB CDF, one can as well in the clustered case, as well as lesser ones in the Weibull and even uniform distributions. As the data become sparser (and gaps become longer and realization scatter increases in the active windows), it is harder to assess visually whether there is actual PWB. Figure 2 illustrates this as we reduce the number of events²²2Due to the random realization, in the PWB case while there are a mean of $N_{\rm total}/4$ events in each activity window, this number varies and so does the total. The three cases for PWB actually have 127, 63, and 30 events. The Weibull distribution case also has 125 rather than 128 events; since the Weibull distribution is well separated from the others, and to enhance clarity, we do not show it for the $N_{\rm total}=64$ and 32 plots. to $N=64$ and then 32.

Since the eye is suspect in determining PWB, we aim to quantify the difference in distributions more rigorously. One common approach to distinguishing between distributions is the relative entropy, or Kullback-Leibler (KL) divergence [19]. For two distributions $p(x)$ and $q(x)$ this is defined as

Generally one does not compare two realized distributions, since they may not have events at the same times $x_{i}$, but rather a realized distribution $p(x_{i})$ is the data and $q(x)$ is the test or model distribution.

We can now start to see why KL divergence is not ideal for our purposes: we cannot take $q(x)$, i.e. the model, to be the PWB distribution, since the important presence of no-event gaps gives $q(x)=0$. There are ways of getting around this using the CDF $Q(x)=\int dx\,q(x)$ rather than the probability density $q(x)$ itself [20], but the results are not wholly satisfactory. Briefly, we do find the ability to distinguish each of the distributions from the uniform random case, and from each other. However, the quantification of the degree of difference is not easily interpreted, and the information of the gaps – the dog not barking – is diluted. (Note that when $p(x_{i})=0$ there is no contribution to the KL divergence, regardless of the model distribution $q(x)$).

Therefore, we adapt a method used in galaxy spatial clustering, where the cosmic web of structure – both connectivity and void regions – carries important information. The friends of friends (FOF) method [21, 22] defines clusters of activity where neighbors lie within a linking length $b$. For our one dimensional time series, this is trivial to implement: $t_{i}-t_{i-1}\leq b$. This is fast and easy to apply to all our distributions. We set the linking length several times larger than the average uniform separation so that clusters will not be (rarely at least) falsely identified in uniform distributions. For example, $b=5/\bar{n}$, where $\bar{n}$ is the number of events divided by the time interval, i.e. the mean density or reciprocal of the inter-event spacing for a uniform distribution, roughly corresponds to a signal to noise $S/N\approx 5$ distinction from a uniform distribution. Empirically, we find this works well.

Friends of friends, like any clustering method, will have difficulties if there are few data points, but note that even if there is only one event in a window, FOF will recognize it as its own cluster unless $\bar{n}$ is too small. In any case, we will end up using FOF as a means of identifying windows, and turn in Sec. III to other statistical techniques for robust quantification of their characteristics. In the end, we will employ FOF quantitatively simply as a useful guide to reasonable priors for a detailed estimation procedure.

An immediate output of applying the FOF method are the values from each activity cluster found of the length of activity windows $a_{j}$, of gaps $g_{j}$ between them, and the possible periods from summing consecutive active and inactive times $T_{j}=a_{j}+g_{j}$ and the active fractions or duty cycles $f_{j}=a_{j}/T_{j}$. For true PWB we expect consistency (i.e. a narrow distribution) in $a_{j}$, $g_{j}$, etc. while these would be widely scattered or sparse for event distributions without some periodic behavior.

Since PWB events occur somewhere in an active window, not necessarily at the extremes, we expect the measured values of $a_{j}$ to give lower limits on the true value $a$, measured values of $g_{j}$ to give upper limits on the true $g$, and the period $a_{j}+g_{j}$ to have some scatter around the true period $T$. All estimates become more accurate with more events in the time series; Section III discusses accurate characterization of the PWB.

Figure 3 shows the results of the FOF clustering analysis on the four types of distributions realized. The uniform random distribution was found to have two, disparate clusters, one extending from the first event to 34% of the observing duration, the other extending for the last 60%, showing no clear periodic windowing. The Weibull distribution has four clusters, with activity window lengths scattered from below 1% to 64% of the duration, again no PWB evident. The clustered distribution has five clusters with three activity windows at 2% length, but others at 7% and 17% length; the related “periods” (sums of consecutive active and gap states) are more diverse, from 8% to 33% length, so no PWB is falsely detected. To compare the different distributions more clearly, we normalize each by their maximum values of period and activity found, $T_{\rm max}$ and $a_{\rm max}$ respectively.

Finally, analysis of the actual PWB distribution case has three activity windows of 15% length and one of 11% length, with periods grouped from 23% to 27%. Recall that the input for generating the random realization was activity 0.15 and period 0.25, so FOF successfully reconstructs the truth. One can quantify this by computing the mean and standard deviation of the activity window length, for example, to assess the peaked nature of the estimation distribution³³3We use a weighted mean, $$\bar{a}=\frac{\sum_{j}N_{j}a_{j}}{\sum_{j}N_{j}}\,,$$ (2) where $N_{j}$ is the number of events within the activity window of duration $a_{j}$. To the extent the scatter goes as $1/\sqrt{N_{j}}$, this is inverse variance weighting. Using this mean, we calculate the standard deviation..

We go into characterization of PWB properties in more detail in the next section, going beyond FOF. Here we have motivated that PWB can be recognized, and that other distributions, even those that have innate clustering (Weibull or the correlation function clustered cases) do not lead to false PWB identification.

Now that we have some confidence that PWB can be identified accurately, we turn to robust characterization of PWB properties such as the period and activity fraction.

Real observational data can be more difficult, and scientifically rewarding. We apply our methods to actual measurements of the source SGR1935+2154, testing whether it exhibits PWB, and giving robust estimates of its period and active window fraction. These characteristics can constrain models of the origin of the emission and properties of the system.

SGR1935+2154, a magnetar within our Galaxy, was identified as the source of two FRBs occurring UT 2020 April 28 [9, 10]. The local nature makes measurements of the object and its environment much easier and with much greater detail. In fact, a linkage between the production of soft gamma bursts in this source and FRBs is naturally suggested: the same two events as the FRBs, with the appropriate delay due to dispersion by interstellar electrons, were detected by $\gamma$-ray instruments [23, 24, 25, 26, 27] (though some authors suggest that different source types may be responsible for these vs extragalactic FRBs [28]).

As an example of the constraining power of these measurements, in [8] the SGR1935+2154 soft gamma burst PWB period was found to be 231 days, but the binary comb model for FRB [29] is limited to periods $\lesssim 100$ days [30]. Such a model may perhaps be eliminated for SGR1935+2154 or other PWB long-period magnetars [8], pointing to other models such as isolated neutron star precession (e.g. [31]) as the source of periodicity, coupled with a nonperiodic emission mechanism. More generally, given an excellent knowledge of the PWB parameters, correlating additional observational properties with these parameters may shed further light on the physical mechanisms at play. For example, with robust knowledge of the phase of the window boundaries, one could look not only for correlations with the modulation period, but also for effects at the window boundaries in SGR burst intensity, fluence, and spectral characteristics. Comparison of the variation in soft gamma burst polarization measurements (coming from the next generation of instruments) with the PWB period and the window boundaries could shed light on the role of the magnetic field in these modulations, as well as in the role of magnetic field orientation in the emission mechanism.

Long term monitoring of SGR1935+2154 has provided several years of data from a number of $\gamma$-ray instruments. In [8], for some range of periods the data were folded at various trial periods, and the activity fraction – the fraction of a period that would be consistent with all event data, was calculated, and the period with the minimum activity fraction was taken to be the “best” period. That analysis did not provide an examination of the uncertainty in the period or active fraction. Its conclusion that there was actually a periodic windowed behavior arose from comparison to uniform random events. Here we use different methods and quantification, and a clustering analysis.

We apply our methods to burst data from the Third Interplanetary Network (IPN3 [32]). IPN3 includes numerous spacecraft with X-ray and gamma-ray sensitive instruments, but notably the Konus instrument on the Wind spacecraft. This instrument is in orbit around the sun at Lagrange point 1, far from earth, and so provides a nearly continuous, unobstructed view of the entire sky, and a more constant background than for low-earth orbit instruments. These are ideal properties for time series monitoring⁴⁴4While there is some heterogeneity in instruments and coverage, the main instruments other than Konus are in low earth orbit, without any long-term changing viewing zones and hence no bias for or against periods in the hundreds of days range. An analysis for particular sets of instruments was carried out in [8] and found results consistent with each other and the results here. .

Carrying out our procedure for identifying and characterizing PWB, we first examine the CDF, as in Fig. 1, but here applied to SGR1935+2154 data (159 events identified in the IPN3 SGR list for SGR1935+2154 as of 2021 February 1 [32]). Figure 6 shows the resulting CDF. While there are clearly episodes of activity and gaps of inactivity, it is difficult to tell by eye if there is PWB.

Applying the FOF method, we carry out the period analysis in Figure 7. While the activity window lengths and gaps scatter greatly, the pseudoperiods (sum of consecutive active window lengths and gap lengths) show an interesting pattern. We exhibit $T/T_{\rm max}$, and see a concentration around $T/T_{\rm max}\sim 0.2$, or $T\sim 0.1$.

The FOF method gives the estimate for the period (not pseudoperiod ratio) of $T=[0.056,0.129]$. While a broad estimation, due to some empty activity windows increasing the uncertainty, it still provides a useful prior for the more incisive likelihood analysis. In addition we find the activity window length $a=[0.049,0.077]$ and six active windows.

We proceed with our likelihood analysis and parameter determination of the SGR1935+2154 data in the same manner as the simulated cases of Sec. III. The final results determine a global optimal period $T=0.1074$ and active fraction $a/T=0.554$ (i.e. $a=0.0595$). Converting back to days by rescaling to the duration of observations, this implies the PWB has period $T=230.6$ days, with an active fraction of 55.4%. Figure 8 shows the results, along with the local optima boundary. Using the efficiency criterion, analysis of the SGR1935+2154 data yields a range of $T\in[0.1066,0.1109]$, or $[-0.7\%,+3.2\%]$ uncertainty on the period.

Figure 9 presents the estimation of the phase parameter, with a best estimate of $\tau/a=-0.101$, i.e. the first detected burst of the time series occurred 10% of the way through the activity window. The estimation of this is very tight ($-0.2\%,+0.9\%$) at $(a/T)_{\rm min}$, broadening as shown in the figure as one moves away from the local optima (here we show a shift by 0.02 from the minimum $a/T$, the same fractional difference as the 0.01 shift used in Fig. 5).

Finally, we note that the FOF estimations for the period $T$ and active window length $a$ do include the likelihood optimization results of $T=0.1074$, $a/T=0.554$ (i.e. $a=0.0595)$, and the optimization results lie well inside the prior information from the FOF analysis. Our results agree as well with those from [8].

Astrophysical bursts occur in observations throughout the electromagnetic spectrum, from the radio to optical to gamma ray. Repeated outbursts indicate the source is not totally disrupted, and periodic bursting points at some physics connected with, e.g., rotation or an orbital companion. An intriguing middle ground that is becoming more recognized with further data is periodic windowed behavior, where activity windows, rather than the burst events themselves, have periodicity. This can also provide important clues to the astrophysical mechanism of the burst, and system characteristics.

Analysis methods for strict periodicity often fall short when dealing with PWB, as the duty cycle is important, activity windows can be empty of events, and the time series of events can be distributed in a complicated manner. We have emphasized that the lack of bursts carries critical information that must fold into the analysis, and we develop a time domain method that takes this into account.

The cumulative distribution function of event intervals works well at identifying whether or not PWB is a reasonable possibility. We test this technique for four distributions: uniform random, Weibull, PWB, and a special distribution with correlated clustering. Given reasonable indications identifying PWB from the CDF analysis, we then draw on the friends of friends technique from galaxy clustering to characterize the PWB. This FOF analysis delivers quantitative estimates of the period $T$, active window length $a$ (and hence duty cycle $a/T$), and observing phase $\tau$. For our test cases of mock data, the estimates accurately reconstruct the input and have $\sim 3-5\%$ precision for data with at least 100 events and $\gtrsim 4$ populated windows. However, we view the FOF analysis as a guide toward carrying out a full likelihood analysis (where the FOF can serve in setting reasonable priors).

For the likelihood analysis we use an optimization grid approach, due to the low dimensionality of the parameter space, increased accuracy, and that the posterior surface is actually a broad plateau, not an isolated peak. However we do find consistent results with both a simulated annealing approach and a Markov Chain Monte Carlo approach. We minimize the duty cycle that agrees with the data, and find accurate estimates of the period to $\lesssim 1\%$ uncertainty on the mock data. Analyzing real observational data on the source SGR1935+2154 we first identify that PWB is reasonable for the data, and then characterize it as having period $T=230.6$ days, with $\sim 1.9\%$ uncertainty, and duty cycle $a/T=55.4\%$.

Finding PWB for what is truly a random distribution is highly unlikely: if we consider the probability that a uniform distribution realization would avoid all the 10 spans of time that PWB predicts no activity (let alone have the active windows in a periodic pattern), this is $P=(1-0.554)^{10}\approx 3\times 10^{-4}$. For a truly uniformly random distribution one could take into account not just the number of windows but the total number of events, so the probability would be $P\approx(1-a)^{N}\approx 2^{-159}\ll 1$. Nevertheless, the ultimate proof will be predictivity: if the values for the period and active fraction derived above are correct, the next two active windows (which admittedly are not guaranteed to have activity) are from June 1 - October 7, 2021 and January 18 - May 26, 2022 and we predict no activity outside of our active windows⁵⁵5On June 24, 2021 a burst was detected by the Fermi Gamma-Ray Space Telescope, in our predicted window. Equally importantly, no bursts were detected outside our predicted windows..

Numerous next generation time domain surveys in a wide range of wavelength bands (e.g. LSST [33] in the optical, DSA-2000 [34] and CMB-S4 [35] in the radio and submillimeter, wide-field instruments such as STROBE-X WFM [36] and SVOM ECLAIRs [37] in the X-ray/gamma ray bands) will greatly increase the database and diversity of repeating sources with possible PWB. The efficient methods presented here give a straightforward path for analysis, identification as PWB (vs, e.g., simply clustering), and its characterization. Accurate estimations of PWB, and the period, duty cycle, and phase, offer the potential for significant advances in understanding the physics of energetic bursts and the properties of the repeating outburst systems, in a wide variety of astrophysical contexts.