Comprehensive Analysis of a Dense Sample of FRB 121102 Bursts

Characterizing the diversity and event rates of FRBs as a function of observing frequency is critical for understanding their nature, the extreme emission physics responsible for FRB and pulsar emission, and the relationship between these two classes of objects. Many surveys have sought comprehensive estimates of these values, all using different observing frequencies, telescopes, and search algorithms but without characterizing the completeness of their search. The Petabyte Project ⁶ (TPP) aims to address these issues to provide robust event rate estimates and discoveries in several petabytes of new and archival radio data. TPP will perform a uniform search for FRBs in an unprecedented amount of archival data and better probe transients closer, farther, and at higher radio frequencies than previous searches. Our search will have a robust internal assessment of completeness, allowing us to confidently project the frequency-dependent rates of FRBs and other transients.

TPP will use Your (a recursive acronym for "your unified reader"; Aggarwal et al. 2020) to ingest the data and heimdall (Barsdell 2012) ⁷ to search it for single pulses. The deep learning based classifier fetch (Agarwal et al. 2020a) ⁸ will then used to classify the candidates identified by heimdall. The data will be searched up to a DM of 5000 pc cm⁻³ (or more, if possible) and a pulse width of 32 ms. This pipeline can easily be modified to search for higher DMs and pulse widths, in specific cases. All of the candidates above a signal-to-noise ratio (S/N) of 6 classified as astrophysical by fetch will be manually verified. The maximum DM and the pulse width to be searched is governed by the observing frequency, data resolution, and GPU memory and hence would be dealt with on a case-by-case basis. The data and results presented in this paper were analyzed under TPP, using an early version of the TPP pipeline. In the following subsection, we discuss the details of the single-pulse search pipeline used in this analysis.

Within Your, we used your_heimdall.py, which runs heimdall on the psrfits data files collected by PUPPI for the single-pulse search. We used two search strategies: (1) a DM range between 450–650 pc cm⁻³ with a DM tolerance ⁹ of 1%, and (2) a DM range between 10–5000 pc cm⁻³ with a tolerance of 25%, with a maximum pulse width of 84 ms in both cases (note that the widest FRB 121102 pulse reported at this frequency had a width of ∼35 ms; see Cruces et al. 2020). The searches resulted in 1428 and 11,276 candidates, respectively. For each candidate, we extracted a segment of the data (which we hereafter refer to as a "cutout") centered at the arrival time (referenced to the top of the observing band) as reported by heimdall with a time window equal to twice the dispersion delay using your_candmaker.py. We then used spectral kurtosis RFI mitigation with a 3σ threshold (Nita & Gary 2010) to identify and excise frequency channels corrupted by RFI (see Figure 2 for the fraction of bursts for which a frequency channel was flagged due to RFI) and used this cleaned data to create dedispersed frequency−time images where a factor of width/2 was used to decimate the time axis. We created the DM–time image by dedispersing the data from zero to twice the reported DM and simultaneously decimating the time axis as above (for more details on the candidate pre-processing, see Agarwal et al. 2020a). These cutouts are then used by fetch to label FRBs and RFI, and were also manually verified. In total, we found 133 bursts with DMs consistent with that of FRB 121102 (i.e., between 550 and 580 pc cm⁻³) with 93 new bursts as compared to the previously published results. We highlight some important differences between our single-pulse search pipeline and the one used by Gourdji et al. (2019) in Section 5.6 to explain the new burst detections. We did not detect any bursts at other DMs. Our search missed one burst, B33, reported by Gourdji et al. (2019), probably because it was weak and narrowband (see Section 5.7 for caveats of our search). Figure 1 shows the dynamic spectra of some of the bursts. Candidate cutouts for all of the bursts are available on GitHub. ¹⁰

Zoom In Zoom Out Reset image size

We define the completeness limit as the pulse energy (also known as fluence) value above which any burst emitted during the observation would be detected. Determining the completeness limit of any single-pulse search is, therefore, essential to defining the sample of bursts to be used for statistical analyses. The most robust method of determining the completeness limit involves an exhaustive injection analysis. In such analyses, simulated transients (with varying properties) are injected into background data, and by analyzing the transients that were recovered (or missed), one can determine the completeness limit of a search (Farah et al. 2019; Agarwal et al. 2020b; Gupta et al. 2021; Li et al. 2021; The CHIME/FRB Collaboration et al. 2021). Such an analysis requires access to a large amount of native-resolution data observed with the same telescope and observing configuration as the search data.

Zoom In Zoom Out Reset image size

As we had access to decimated data for just two observations, we could not do such an injection analysis. We, therefore, estimate the completeness limit from the radiometer equation. We use a conservative approach by including the effect of RFI mitigation, which reduces search sensitivity. We flagged more than 35% of data for many candidates due to RFI, leading to a usable bandwidth of 500 MHz. Using this smaller bandwidth and a nominal pulse width of 1 ms, the fluence limit above an S/N of 8 is 0.0216 Jy ms. We use bursts with a fitted fluence above this limit in the burst sample analysis.

To measure the properties of the bursts in our sample, we perform spectro-temporal modeling of all of the detections. We model each component of the bursts' spectra with a Gaussian function and the profile using a Gaussian convolved with a one-sided exponential function to represent the Gaussian pulse and the scattering tail. Therefore, we model each component of the burst (at an observing frequency f and time t), using the following function,

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal F }(f,t;S,{\mu }_{f},{\sigma }_{f},{\mu }_{\mathrm{DM}},{\sigma }_{t},{\tau }_{\mathrm{sc}})\\ \quad =\,S\times 2.355{\sigma }_{f}\\ \quad \times \,{{ \mathcal P }}_{f}(t;{\mu }_{\mathrm{DM}},{\sigma }_{t},{\tau }_{\mathrm{sc}})\\ \quad \times \,{ \mathcal G }(f;{\mu }_{f},{\sigma }_{f}).\end{array}\end{eqnarray} \tag{ 1 }$$

Here, S is the fluence of the component, ${ \mathcal G }$ is the Gaussian function to model the spectra

$$\begin{eqnarray}&&{ \mathcal G }(x;{\mu }_{x},{\sigma }_{x})=\displaystyle \frac{1}{{\sigma }_{x}\sqrt{2\pi }}\exp \left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\left(x-{\mu }_{x}\right)}^{2}}{{\sigma }_{x}^{2}}\right)\end{eqnarray} \tag{ 2 }$$

and ${ \mathcal P }$ is used to model the profile. For bursts with scattering, we use the Gaussian convolved with a one-sided exponential function (McKinnon 2014),

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal P }}_{f}^{\mathrm{scat}}(t;{\mu }_{\mathrm{DM}},{\sigma }_{t},{\tau }_{\mathrm{sc}})\\ \quad =\,\displaystyle \frac{1}{2{\tau }_{\mathrm{sc}}}\left\{1+\mathrm{erf}\left[\displaystyle \frac{t-({\mu }_{\mathrm{DM}}+{\sigma }_{t}^{2}/{\tau }_{\mathrm{sc}})}{{\sigma }_{t}\sqrt{2}}\right]\right\}\\ \qquad \times \,\exp \left(\displaystyle \frac{{\sigma }_{t}^{2}}{2{\tau }_{\mathrm{sc}}^{2}}\right)\exp \left(\displaystyle \frac{t-{\mu }_{\mathrm{DM}}}{{\tau }_{\mathrm{sc}}}\right),\end{array}\end{eqnarray} \tag{ 3 }$$

where the mean of the Gaussian pulse after accounting for the dispersion delay at that respective frequency

$$\begin{eqnarray}&&{\mu }_{\mathrm{DM}}={\mu }_{t}-4.148808\times {10}^{3}\,\mathrm{DM}\left(\displaystyle \frac{1}{{f}^{2}}-\displaystyle \frac{1}{{f}_{\mathrm{top}}^{2}}\right)\,{\rm{s}}.\end{eqnarray} \tag{ 4 }$$

Here, f_top is the highest frequency in the band (in MHz), f is the frequency of a channel (in MHz) and μ_t is the mean of the Gaussian pulse at f_top. At each frequency channel, the scattering timescale

$$\begin{eqnarray}&&{\tau }_{\mathrm{sc}}(f)={\tau }_{\mathrm{sc}}{\left(\displaystyle \frac{f}{{f}_{\mathrm{ref}}}\right)}^{-4}.\end{eqnarray} \tag{ 5 }$$

In this expression, f_ref is the reference frequency and is set to be 1 GHz. The exponent (–4) is assuming a normal distribution of plasma-density inhomogeneities. Finally, we defined ${{ \mathcal P }}_{f}(t;{\mu }_{\mathrm{DM}},{\sigma }_{t},{\tau }_{\mathrm{sc}})$ to be a Gaussian ${ \mathcal G }(t;{\mu }_{\mathrm{DM}},{\sigma }_{t})$ for τ_sc/σ_t < 6, and ${{ \mathcal P }}_{f}^{\mathrm{scat}}(t;{\mu }_{\mathrm{DM}},{\sigma }_{t},{\tau }_{\mathrm{sc}})$ for τ_sc/σ_t > 6. We used this value for the cutoff in order to maintain numerical stability while calculating Equation (3).

It follows from the above discussion that our model ${ \mathcal F }$ is generated using seven parameters: S, μ_f , σ_f , μ_t , σ_t , τ_sc, and DM. Using this model we fit the burst spectrograms, as described in the next section.

While fitting for complex FRB bursts is an arduous task, we scrupulously automate the entire procedure and present it as the Python package burstfit. ¹¹ burstfit provides a framework to model any spectrogram consisting of any complex FRB or pulsar pulse using robust methods. It can easily incorporate any user-defined Python function(s) to model the profile, spectra, and spectrogram and is not limited to the functions we have implemented for this current analysis. burstfit primarily consists of the following five steps. Figure 3 shows the flowchart with the various stages of fitting in BURSTFIT.

3.5.1. Data Preparation

3.5.2. Stage 1: Single-component Fitting

3.5.3. Stage 2: Statistical Tests

3.5.4. Stage 3: Multicomponent Fitting

In cases where multiple components were found, we performed another stage of the fitting. Here, we generated a combined model consisting of all of the components and fit for all components by comparing our model with the original spectrogram. This combined model (${{ \mathcal F }}_{\mathrm{all}}$) was generated by summing together the individual component models (${{ \mathcal F }}_{i}$) for all N components,

$$\begin{eqnarray}&&{{ \mathcal F }}_{\mathrm{all}}=\sum _{i}^{N}{{ \mathcal F }}_{i}.\end{eqnarray} \tag{ 6 }$$

Again, we used scipy.curve_fit for fitting. Here, the fit results of the individual component fits from previous stages were used as the initial guess for the parameters in curve_fit.

3.5.5. Stage 4: MCMC

While scipy.curve_fit is sufficient for fitting in many scenarios, in our testing we found that the estimates and the errors reported by scipy.curve_fit were not robust for our purposes. In many cases, the errors reported by scipy.curve_fit were possibly underestimated, and fitted results were highly susceptible to the choice of input parameter bounds. This was especially true for low-significance bursts and multiple-component bursts, where the least-squares-minimization technique struggles to find a good solution. Therefore, we added another stage to our fitting procedure and used a Markov Chain Monte Carlo (MCMC) to obtain the final fitting results. We used the results of previous stages (that used scipy.curve_fit) as initial estimates to determine the starting positions of the walkers for the MCMC. An advantage of the MCMC procedure is that it provides the full posterior distribution of all of the fitted parameters, which we could then use to estimate the errors and further follow-up analysis of the burst sample. We used the Goodman and Weare affine invariant sampler (Goodman & Weare 2010) as implemented in emcee (Foreman-Mackey et al. 2013). We used uniform priors for all of the parameters, with the ranges of the priors given in Table 1. We used the log-likelihood function

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }=-0.5\sum {\left(\displaystyle \frac{{ \mathcal S }-{{ \mathcal F }}_{\mathrm{all}}}{\sigma }\right)}^{2},\end{eqnarray} \tag{ 7 }$$

where ${ \mathcal S }$ refers to the original spectrogram, ${{ \mathcal F }}_{\mathrm{all}}$ refers to the model, and σ is the off-pulse standard deviation of the measured spectrogram. The sum is over all of the pixels in the two spectrograms. We used an autocorrelation analysis to determine when the MCMC converged. ¹³ We then estimated the burn-in ¹⁴ using the autocorrelation time and used the remaining samples to determine the fitting results. To decide if scattering was present in a burst, we used the percentage of samples with τ_sc/σ_t < 6. If this percentage was greater than 50%, we concluded that scattering was not present (or was very small) in that burst. We do not report scattering timescales for such bursts.

Table 1. Priors Used in the MCMC Fitting

Parameter	Minimum	Maximum
S	0	500 × max(time_series)$\times {\sigma }_{t}^{\mathrm{fit}}$
μf	−2 × Nf	3 × Nf
σf	0	5 × Nf
μt	$0.8\times {\mu }_{t}^{\mathrm{fit}}$	$1.2\times {\mu }_{t}^{\mathrm{fit}}$
σt	0	$1.2\times ({\sigma }_{t}^{\mathrm{fit}}+{\tau }_{\mathrm{sc}}^{\mathrm{fit}})$
τsc	0	$1.2\times ({\sigma }_{t}^{\mathrm{fit}}+{\tau }_{\mathrm{sc}}^{\mathrm{fit}})$
DM	0.8 × DMfit	1.2 × DMfit

Note. Superscript fit refers to the values obtained using fits done in previous stages. N_f refers to the number of frequency channels. The term time_series refers to the 1D array obtained by summing the dedispersed cutout spectrogram along the frequency axis. Subscripts t and f are used for profile and spectra parameters.

Download table as:  ASCIITypeset image

We generated corner plots and fit-result plots to verify the quality of the fits, as shown in Figures 4 and 5. We provide all of the results (output parameters, corner plots, fitting-result plots, etc.) from our analysis in a GitHub repository. ¹⁵

3.5.6. Handling Data Saturation

The data we use in this analysis were recorded as 8-bit unsigned integers. Hence, the data range lies between 0–255, and any signal brighter than 255 is clipped at this value. We noticed data saturation for two bursts (B6 and B121), and hence this effect has been incorporated in our burst modeling. The spectrograms are subtracted by the off-pulse mean (μ_off) and divided by the off-pulse standard deviation (σ_off). While making the spectro-temporal model, we clip the values greater than (255 − μ_off)/σ_off. This effect is visible in Figure 6 for burst B121 where the red dotted–dashed curve and green dotted curve show the fit to the burst spectra with and without clipping, respectively. The fit performed without considering the saturation underestimates the burst's spectral width, leading to an underestimated burst energy.

Zoom In Zoom Out Reset image size

3.5.7. Caveats to Our Fitting Analysis

There are some caveats to our fitting procedure that are worth noting here. First, as with any model-dependent fitting, our analysis and results are dependent on the choice of the functions we use to model the data. We described those functions and our motivation for using them in Section 3.4, but these are not the only proposed methods to model the spectrogram of an FRB. The spectrum of an FRB can also have power-law-like behavior (CHIME/FRB Collaboration et al. 2020), and some other pulse broadening functions can also be used besides the exponential tail used to model the scattering effect (Bhat et al. 2004 ).

Second, the emission of some bursts in our sample was present only in the top part of the band. In many such cases, the emission appeared similar to the tail of a Gaussian function, with the mean lying outside our observing band. Although we allowed for such a mean value to be estimated with the MCMC procedure, our fitting results for such bursts would inevitably be unconstrained. We therefore mark such bursts with a † in Table 2 to highlight this.

Table 2. Properties of the First 10 Bursts

Burst a	μf	σf	S	μt b	σt	τc	DM
ID	(MHz)	(MHz)	(Jy ms)	(MJD)	(ms)	(ms)	(pc cm−3)
B1*	${1560}_{-30}^{+30}$	${210}_{-40}^{+40}$	${0.09}_{-0.02}^{+0.02}$	57,644.408906976(1)	${0.0}_{-0.02}^{+0.02}$	${1.9}_{-0.3}^{+0.3}$	${565.3}_{-0.4}^{+0.4}$
B2*	${1200}_{-10}^{+10}$	${50}_{-10}^{+10}$	${0.043}_{-0.007}^{+0.007}$	57,644.40956768(1)	${1.35}_{-0.05}^{+0.05}$	-	${562.4}_{-0.8}^{+0.8}$
B3.1 †	${2900}_{-600}^{+300}$	${800}_{-200}^{+300}$	${0.6}_{-0.3}^{+0.7}$	57,644.409673699(3)	${0.4}_{-0.2}^{+0.2}$	${1.3}_{-0.7}^{+0.7}$	${566.8}_{-0.9}^{+0.8}$
B3.2 †	${1100}_{-1400}^{+300}$	${1000}_{-1000}^{+2000}$	${0.09}_{-0.04}^{+0.13}$	57,644.40967384(2)	${0.3}_{-0.2}^{+1.4}$	${0.3}_{-0.1}^{+0.7}$	${564.7}_{-0.4}^{+1.5}$
B4 †	${3100}_{-600}^{+200}$	${550}_{-110}^{+80}$	${2}_{-2}^{+4}$	57,644.410072889(4)	${1.1}_{-0.2}^{+0.2}$	${0.3}_{-0.1}^{+0.2}$	${564}_{-1}^{+1}$
B5 †	${2100}_{-1600}^{+900}$	${2700}_{-1200}^{+900}$	${0.19}_{-0.04}^{+0.05}$	57,644.410157834(4)	${0.7}_{-0.3}^{+0.3}$	${1.0}_{-0.3}^{+0.4}$	${562.1}_{-0.6}^{+0.7}$
B6.1	${1393}_{-7}^{+7}$	${183}_{-7}^{+7}$	${0.47}_{-0.03}^{+0.02}$	57,644.411071954(1)	${1.09}_{-0.04}^{+0.03}$	-	${562.3}_{-0.1}^{+0.2}$
B6.2	${1417}_{-5}^{+4}$	${102}_{-4}^{+5}$	${0.33}_{-0.02}^{+0.03}$	57,644.4110719755(9)	${0.57}_{-0.02}^{+0.03}$	-	${560.9}_{-0.1}^{+0.2}$
B7.1 †	${3100}_{-300}^{+200}$	${430}_{-80}^{+60}$	${10}_{-10}^{+20}$	57,644.412240214(5)	${0.7}_{-0.4}^{+0.3}$	${0.8}_{-0.4}^{+0.4}$	${569}_{-3}^{+3}$
B7.2 †	${1460}_{-20}^{+20}$	${90}_{-20}^{+30}$	${0.09}_{-0.01}^{+0.01}$	57,644.41224043(2)	${1.9}_{-1.0}^{+0.9}$	${1.2}_{-0.5}^{+0.6}$	${569}_{-2}^{+3}$
B8 †	${3000}_{-600}^{+200}$	${700}_{-100}^{+100}$	${1.1}_{-0.7}^{+1.5}$	57,644.414123628(4)	${1.0}_{-0.3}^{+0.3}$	${0.8}_{-0.3}^{+0.4}$	${567.5}_{-0.7}^{+0.9}$
B9*	${1430}_{-10}^{+10}$	${75}_{-9}^{+9}$	${0.076}_{-0.003}^{+0.003}$	57,644.41447161(2)	${2.0}_{-0.6}^{+0.6}$	${0.4}_{-0.6}^{+0.6}$	${564}_{-2}^{+2}$
B10	${1630}_{-10}^{+10}$	${82}_{-8}^{+8}$	${0.1}_{-0.01}^{+0.01}$	57,644.414475391(7)	${1.4}_{-0.3}^{+0.3}$	${0.5}_{-0.2}^{+0.3}$	${562}_{-3}^{+3}$

Notes. See Table 5 in the Appendix for the full table. The 1σ errors on the fits are shown in the superscript and subscript of each value in the table. For μ_t , the error on the last significant digit is shown in parenthesis.

^a Burst IDs are chronological. Individual component number (N) for multicomponent bursts are appended to the burst IDs. Bursts modeled only using curve_fit are marked with *. Note that the errors on these bursts could be unreliable and may be either under- or overestimated. Bursts that extend beyond the observable bandwidth can also have unreliable estimates of spectra parameters and fluence (see Section 3.5.7). We mark those bursts with † to indicate that their fluence and spectra parameters could be unconstrained. ^b The term μ_t is the mean of the pulse profile in units of MJD. This can be considered as the arrival time of the pulse. It is referenced to the solar system barycenter, after correcting to infinite frequency using a DM of 560.5 pc cm⁻³. ^c Here, τ is referred to 1 GHz.

Download table as:  ASCIITypeset image

Third, in some cases, the MCMC procedure was unable to find a robust solution. This was due to the presence of RFI, which could not be removed as the FRB signal coincided with the channels heavily corrupted by RFI. It also occurred when there was significant baseline variation in the data close to a weak FRB pulse. This could dominate the MCMC likelihood estimate, and therefore the procedure could not converge on a solution. For such cases, we only used scipy.curve_fit to perform the fits, and we modified the fitting bounds to obtain a visually good fit. We highlight these in Table 2 with *, and note that the values could be unreliable.

We used the converged sample chains from the MCMC fitting for each burst to generate a cumulative corner plot with properties of the whole burst sample. To do this, we randomly selected 1000 samples from the final 25% of the MCMC chains and then concatenated such samples from all of the bursts. We then generated a corner plot using these samples, as shown in Figure 7. ¹⁶ We can now use Figure 7 to infer trends in various spectro-temporal properties of FRB 121102. Table 3 shows the summary statistics of all of the bursts obtained from this analysis. We can see that the spectra of the bursts typically peak around 1650 MHz, and there is a dearth of burst emission below 1300 MHz. Most of the burst spectra peak within the top part of our observing band (i.e 1550–1780 MHz). This behavior possibly extends further to higher frequencies, as is evident from many burst spectra that increase toward the top part of the observing band with their spectral peak possibly lying outside our observing band. Interestingly, Platts et al. (2021) have recently reported complex bifurcating structures in some FRB 121102 bursts below 1250 MHz using higher-resolution data. It is, therefore, possible that the emission of FRB 121102 shows a different behavior below these frequencies, which might also vary with time.

Zoom In Zoom Out Reset image size

Table 3. Results from the Burst Sample Analysis

Parameter	Units	Value
μf	MHz	${1608}_{-200}^{+100}$
σf	MHz	${102}_{-30}^{+130}$
S	Jy ms	${0.13}_{-0.05}^{+0.13}$
σt	ms	${1.1}_{-0.5}^{+0.9}$
τ	ms	${0.7}_{-0.4}^{+0.9}$
DM	pc cm−3	${564}_{-3}^{+5}$

Note. The values represent the median values with 1σ errors.

Download table as:  ASCIITypeset image

As already noted by Gourdji et al. (2019), FRB 121102 shows a variety of spectral widths. Using our modeling, we observe that most bursts are frequency modulated and have a typical frequency width of ∼230 MHz. The bursts also show a wide range of intrinsic pulse widths, from 0.4 to 20 ms, and various scattering timescales, up to 3 ms. The median DM of the bursts we observe was 564 pc cm⁻³, with a 1σ variation of ∼4 pc cm⁻³. This variation in DM is also apparent in Figures 7 and 8. This value is consistent with the other published estimates (Gourdji et al. 2019; Cruces et al. 2020; Li et al. 2021; Platts et al. 2021). We also did not see any strong correlation between any two burst properties from Figure 7. Several of the bursts from this sample also show the characteristic subburst drift in frequency during the burst duration, sometimes referred to as the "sad-trombone" effect (Hessels et al. 2019). We did not detect any evidence of upward drifting as predicted by some FRB models (Cordes et al. 2017), and reported by Platts et al. (2021).

Zoom In Zoom Out Reset image size

Many bursts in our sample show multiple components, and we estimated the properties of these components using our fitting procedure. Nine bursts in our sample show two components, while there is one burst with three components (see Table 2). We also note that it is difficult to differentiate between multiple closely spaced bursts and different components from single bursts. This is further complicated by the detection of a very wide (∼35 ms) burst reported by Cruces et al. (2020). As there is no clear consensus on how to resolve this, we visually identified some bursts as components of a nearby burst and reported them as such in Table 2. We consider all of the components as individual bursts for all of the following analysis except the cumulative energy distribution analysis.

Figure 8 shows the scatter plots of various burst properties with respect to the burst time. The bursts from two observations are shown in different colors, and the time is referenced to the first burst of the respective observation. The burst properties do not show any temporal evolution on the seconds-to-minutes timescale. We also did not observe any distinction between the distribution of properties of bursts detected on two consecutive days.

Energy distributions can provide useful intuition into the emission mechanism of the source. Regular pulsar emission typically shows a log-normal distribution, whereas giant pulses show power-law cumulative distributions (Burke-Spolaor et al. 2012). Crab giant pulses have shown evidence that the index depends on pulse width and energy, with flatter indices for weaker and shorter pulses (Popov & Stappers 2007; Karuppusamy et al. 2010; Bera & Chengalur 2019). High-energy magnetar emission has been described by power-law distributions with γ ranging from roughly −1.6 to −1.8 (Cheng et al. 2020). Previous studies of FRB 121102 energy distributions have used a single power-law fit (N(>E) ∝ E^γ ) to model the cumulative distribution and have obtained different values of γ including −0.7 in Law et al. (2017), −1.1 in Cruces et al. (2020), −1.7 in Oostrum et al. (2020), and −1.8 in Gourdji et al. (2019). Another well-studied repeating source, FRB 180916, shows γ = −1.3 at 400 MHz, although recent observations have reported a flattening of the power law at lower energies (CHIME/FRB Collaboration et al. 2020; Pastor-Marazuela et al. 2021).

We calculate the isotropic energy (E) of a burst as

$$\begin{eqnarray}&&\begin{array}{l}E=4\pi {D}_{L}^{2}\times S\,(\mathrm{Jy}\,{\rm{s}})\times 2.355\,{\sigma }_{f}(\mathrm{Hz})\\ \qquad \times {10}^{-23}({\mathrm{erg}}^{-1}{\,\mathrm{cm}}^{-2}{\,\mathrm{Hz}}^{-1}).\end{array}\end{eqnarray} \tag{ 8 }$$

Here, D_L is the luminosity distance to FRB 121102, 972 Mpc, as reported by Tendulkar et al. (2017). S and 2.355 σ_f are the fitted fluence and FWHM of the Gaussian spectra.

To make the cumulative energy distribution, we choose only bursts for which the ±1σ bounds on the spectral peak fell within our observing band. This was done as our fluence estimates obtained from fitting are reliable for bursts within our band and because we are incomplete to the population of bursts that are partially outside our band (Aggarwal 2021). Therefore, from a total of 133 bursts, we obtained 60 bursts that satisfied this criteria. For each of the 60 such bursts, we used the posterior distribution of bandwidths and fluences from the MCMC based fitting analysis to calculate the distribution of energies (using Equation (8)). We then randomly sampled one energy from the burst energy distributions of each of the 60 bursts and generated a cumulative energy distribution using those 60 energies. We repeated this process 1000 times and thereby generated 1000 cumulative energy distributions.

The previous studies of the cumulative energy distribution of FRB 121102 have reported a break in the power law with a flattening toward low energies (Gourdji et al. 2019; Cruces et al. 2020). Also, it was visually evident in our data that the cumulative distribution flattened toward low energies. Therefore, a single power law would not have been sufficient to accommodate the burst energy distribution. Therefore, we used scipy.curve_fit to fit each of these 1000 energy distributions with a broken power law of the form

$$\begin{eqnarray}N(\geqslant E)=\left\{\begin{array}{ll}{E}_{\mathrm{scale}}{\left(\displaystyle \frac{E}{{E}_{\mathrm{break}}}\right)}^{\alpha }, & \mathrm{if}E\lt {E}_{\mathrm{break}}\\ {E}_{\mathrm{scale}}{\left(\displaystyle \frac{E}{{E}_{\mathrm{break}}}\right)}^{\beta }, & \mathrm{if}E\gt {E}_{\mathrm{break}}.\end{array}\right.\end{eqnarray} \tag{ 9 }$$

Here, α and β are the two power-law indices, E_break is the break energy, and E_scale is the energy scaling. Figure 9 shows the cumulative energy distributions above an estimated completeness of 5.8 × 10³⁶ erg, calculated from the aforementioned completeness limit on the fluence, 0.0216 Jy ms (Section 3.3), and a median bandwidth (2.355 σ_f ) of the bursts i.e., 240 MHz. It also shows (in red) the power-law fit to each cumulative energy distribution.

Zoom In Zoom Out Reset image size

The median of the distribution of fitted power-law indices and break energy (with 1σ errors) are given by α = −0.4 ± 0.1, β = −1.8 ± 0.2 and E_break = (2.28 ± 0.19) × 10³⁷ erg. This break at E_break could indicate the actual completeness energy limit of our observations, and we might therefore be incomplete to the bursts with energies below E_break. This could be due to the incompleteness of our observations to the weak, band-limited bursts. The bursts above this energy are well fitted by the power law of index β = −1.8 ± 0.2. The break in energy distribution has also been reported for other repeating FRBs (Pastor-Marazuela et al. 2021). It is worth noting that our higher-energy power-law index β is also consistent with the power-law index estimated by Gourdji et al. (2019) above a completeness threshold of 2 × 10³⁷ erg. In the context of pulse-energy distributions, the similarity of the power-law indices of FRB 121102 with both those of Crab giant pulses and magnetar pulses might also imply a common origin (Lyu et al. 2021).

5.1.1. Testing for a High-energy Break

The left panel in Figure 10 shows the distribution of wait times between bursts, which follows a bimodal distribution as also seen in previous studies (Gourdji et al. 2019; Li et al. 2019). On wait times greater than 1 s, we used scipy.curve_fit to fit a log-normal function to the main distribution, finding a peak at 74.8 ± 0.1 s. The peak of our wait-time distribution is significantly lower than the 207 ± 1 s peak found in these data by Gourdji et al. (2019), due to our increased sample of bursts filling in the wait-time gaps. As more bursts are discovered in an observation of constant length, the average time between bursts decreases, lowering the peak of the wait-time distribution. Our findings most closely match the wait-time distribution of Zhang et al. (2018), which peaks at ∼67 s; while the original paper did not report an exact peak, we used their publicly available data ¹⁷ to perform our wait-time analysis. These similarities suggest that careful single-pulse searches using machine-learning algorithms allow us to obtain a robust sample of bursts that accurately reflects the burst population.

Zoom In Zoom Out Reset image size

As in previous studies (e.g., Katz 2018; Gourdji et al. 2019; Li et al. 2019), we find a smaller population of bursts with subsecond separations. However, unlike the previously reported distributions, which cluster around tens of milliseconds, our subsecond burst separations span the range of tens of milliseconds up to nearly one second without as clear a break between the two distributions. For this analysis, we assumed that each closely spaced pair of bursts is composed of two separate bursts instead of components of a single broader burst, leading to the larger distribution of short wait times compared to other papers. This assumption will not drastically alter the fitted log-normal distribution since the subsecond wait-time population is small, and their removal will not significantly alter the shape of the main distribution. As in Gourdji et al. (2019) and Li et al. (2019), we found that the wait time between bursts and their relative fluences are not correlated.

To quantify the change in the wait-time peak as a function of the number of bursts in a sample, we performed the same wait-time analysis for a random selection of our bursts over a range of burst numbers, as seen in the right-hand panel of Figure 10. Each point represents a fit to 300 random selections of that number of bursts, with the error bars representing the standard deviation of all of the fitted peaks. Our findings show the expected effect that as more bursts are included in a sample from a constant-length observation, the average time between bursts will decrease along with the fitted wait-time distribution peak. We find that the distribution of fitted wait times peaks exponentially decays with added bursts with a timescale of ∼29 s. We observe the same effect in the Zhang et al. (2018) data set, which shows a timescale of ∼25.5 s. The peak obtained by Gourdji et al. (2019) using these data is shown in the figure and matches with our fitted exponential curve. We performed a similar analysis by filtering out the lowest fluence bursts from our sample and found that the wait-time peak increases as the minimum fluence limit is increased and weaker bursts are excluded. This serves to explain the higher wait-time peaks calculated by previous papers with fewer bursts and higher fluence limits.

In this section we discuss the results of various periodicity searches tested on this burst sample to search for a short-period periodicity (millisecond to hundreds of seconds).

5.3.1. Difference Search

5.3.2. Fast Folding Algorithm

5.3.3.  frbpa

5.3.4. Lomb−Scargle

Previous observations of FRB 121102 have found significant evidence for pulse clustering on short timescales, where the burst separations deviate from a Poissonian distribution (Oppermann et al. 2018; Cruces et al. 2020; Oostrum et al. 2020). The Weibull distribution, as described in Oppermann et al. (2018), is a modification of the Poisson distribution, with a shape parameter k describing the degree of pulse clustering. Clustering is present for k < 1 with lower values corresponding to more clustering, while k = 1 reduces to the Poissonian case; a value of k ≫ 1 causes the distribution to peak more sharply at the event rate and indicates a periodic signal. A better understanding of the burst statistics may help us understand the progenitor of FRB 121102 and help strategize the timing of future observations.

Figure 11 shows the cumulative probability density of the wait times between consecutive bursts, fitted to the Weibull and Poisson cumulative density functions as defined by Oppermann et al. (2018). The fitted values are given in Table 4, as well as the reduced chi-squared statistic and the coefficient of determination, r², which ranges from zero to one with a value of one representing a perfect fit. In Figure 11, we also plot the values of the reduced chi-squared statistic for the Poisson and Weibull distributions when fitted to only the wait times longer than a range of chosen minimum wait times. We observe that both the Poisson and Weibull distributions fit the main population of longer wait times much better than the entire set of bursts. The Weibull distribution's ability to account for clustering allows it to have a significantly better fit when including shorter wait times, as its reduced chi-squared statistic is a factor of 18 smaller than that of the corresponding Poisson distribution (see Table 4). However, the Weibull distribution is only slightly favored over the Poisson distribution when fitting only longer wait times. We find that the reduced chi-squared statistic is equal to 0.368 for the Weibull fit to wait times greater than one second, indicating that the distribution is overfitted. These findings may indicate that the main distribution of bursts with longer wait times roughly follows a Poissonian distribution, while the entire burst rate distribution cannot be accurately described with solely a Poisson or Weibull distribution. This may result from our decision to consider each burst as a separate burst rather than a subcomponent of a broad burst.

Zoom In Zoom Out Reset image size

In addition to their observations, Cruces et al. (2020) used the original data set from Gourdji et al. (2019) to study the burst rate statistics and found that the addition or removal of the subsecond wait-time population in each data set significantly impacts the extrapolated burst rate behavior. In each case, removing these short wait times led to the fitted Weibull shape parameter k increasing toward one, further indicating that the main distribution of pulses may follow a Poissonian distribution while the shorter distributions do not. However, the sample used in Cruces et al. (2020) only had two subsecond wait times; our more extensive sample of short wait times allows us to confirm this behavior with more statistical significance. In Figure 12, we plot the fitted burst rates for the Weibull and Poisson distributions and the Weibull shape parameter k as a function of minimum wait time. As the minimum wait time increases, both burst rates converge to a rate of roughly 45 bursts per hour. We also find that the fitted value of k increases with the minimum wait time, reaching a value of k = 1 at a minimum wait-time cutoff of roughly 0.1 s.

Zoom In Zoom Out Reset image size

Based on our results, any progenitor model proposed for FRB 121102 would have to explain the following observations: (1) band-limited emission; (2) varying peak emission of the spectra and its absence below 1300MHz; (3) median scattering timescale of 0.7 ms, with a maximum value of around 2 ms; (4) rapid variability of these three properties at second timescales. Further, some of these observations have also been reported for other FRBs (Shannon et al. 2018; Kumar et al. 2021; Pastor-Marazuela et al. 2021).

In this work, we have presented 93 additional bursts detected upon reprocessing of the data presented in Gourdji et al. (2019). Figure 13 shows the distribution of properties of new bursts detected with our pipeline as compared to the ones already published by Gourdji et al. (2019). We performed Kolmogorov–Smirnov (KS) tests to compare the distributions of fitted parameters of old versus the new bursts. The distribution of S, σ_t , and DM are similar, while those for μ_f , σ_f , and τ are different. This indicates that the searches carried out by Gourdji et al. (2019) missed bursts that span the entire range of these parameter values, rather than just the weaker bursts.

Zoom In Zoom Out Reset image size

Recently, Li et al. (2021) reported a large sample of bursts from FRB 121102 detected using the Five-hundred-meter Aperture Spherical radio Telescope. The mean of the wait-time distribution estimated from that larger sample of bursts (70 ± 12s) is consistent with what we report in Section 5.2. They also did not detect any short-term periodicity, similar to our findings. Further, they report a bimodal energy distribution for FRB 121102. Aggarwal (2021) highlighted that this bimodality disappears when the burst bandwidth, instead of center frequency of the observing band, is used to calculate the energy. Moreover, the burst bandwidths reported in Li et al. (2021) were estimated visually and not using a fitting procedure. This can lead to observational biases that will make the interpretation of intrinsic energy distributions difficult (Aggarwal 2021).

All of the tools and software used in our pipeline are independent of the ones used in the original work (Gourdji et al. 2019). This brings the critical question of understanding why our pipeline detected more bursts or why the original work missed the bursts. While an exhaustive comparison of the two pipelines using a standard data set consisting of simulated FRBs is beyond the scope of this work, here, we try to investigate the reasons for different results based on our understanding of the software used. We discuss three reasons that might contribute to very different recovery rates ¹⁸ of the two pipelines.

5.6.1. Threshold Signal-to-noise Ratio

Although Gourdji et al. (2019) used an S/N threshold of 6 for the search, they discarded candidate groups with a maximum S/N less than 8. We used an S/N threshold of 6 in our search. Therefore, they would have missed low S/N candidates.

Assuming that the higher-energy power-law slope of −1.8 (estimated in Section 5.1) is intrinsic to FRB 121102, we can estimate the expected ratio of the number of bursts with an S/N greater than 6 to that above an S/N of 8. This is given by ${\left(6/8\right)}^{-1.8}=1.7$. The observed ratio of the number of bursts above an S/N of 6 (N = 133) to that above an S/N of 8 (N = 70) is 133/70 = 1.9. This implies that we detected more bursts between S/Ns of 6 and 8 than expected from the power-law distribution of energies. But, it is important to note that this simple estimate relies on the following assumptions:

• 1.  
  Our observations are complete to bursts with an S/N < 8. But the flatter energy distribution at low energies indicates that this might not be true.
• 2.  
  The fluence and energy distributions are similar. This is true only when there is no incompleteness due to banded nature of the burst spectra (Aggarwal 2021).
• 3.  
  The burst energy distribution can be modeled by a single power law, even at low energies.

Further, we now estimate the number of purely noise candidates we expect above an S/N greater than 6 in our search. The chance probability of a purely noise candidate with an S/N greater than 6 is P = 9.9 × 10⁻¹⁰. The number of trials in our search can be calculated by:

$$\begin{eqnarray}&&{N}_{{\rm{t}}{\rm{r}}{\rm{i}}{\rm{a}}{\rm{l}}}={N}_{{\rm{D}}{\rm{M}}}\displaystyle \sum _{i=0}^{10}\displaystyle \frac{{N}_{{\rm{t}}{\rm{i}}{\rm{m}}{\rm{e}}}}{{2}^{i}}.\end{eqnarray} \tag{ 10 }$$

Here, N_trial is the total number of trials in our search, N_DM is the number of trial DMs (65, between 450 and 650 pc cm⁻³ at a tolerance of 1%), and N_time is the total number of time samples in the data (∼3 hr at a time resolution of 81.92 μs). The sum is over the boxcar widths searched (2⁰–2¹⁰, doubling at each step). Therefore, we expect P × N_trial ∼ 20 purely noise events, with an S/N greater than 6, for this number of trials. But the number of events we detected with an S/N > 6 in our single-pulse search was 1427 (with 133 FRB and 1294 RFI candidates).

This shows that we detected many more candidates above an S/N of 6 than expected from pure Gaussian noise, implying that data is non-Gaussian. Most of the candidates we obtained were due to RFI, which is expected. It is possible that some weak events are still misidentified; however, this will not influence the results of our analysis.

5.6.2. Single-pulse Search Software

5.6.3. RFI Mitigation and Classification

5.6.4. General Comments

Finally, it is appropriate to highlight and reiterate four main caveats to the analysis presented in this work. (1) The data used in this analysis were downsampled to 64 frequency channels. We did not have access to the native-resolution data; therefore, all of the search and analysis was performed on downsampled data. The sensitivity of a single-pulse search would be higher on native-resolution data; therefore, our pipeline may have missed some pulses. (2) We only performed a search on the data that averaged over the full bandwidth to create time series. Given the band-limited nature of many bursts, a subband search on these data might reveal weaker and narrower bursts. (3) As mentioned previously, our estimate of the completeness limit is not robust, as such a robust estimate requires an injection analysis that was not possible with the available data. (4) The reported properties of the bursts depend on the assumption that the burst spectra and profile follow the assumed functional forms used for fitting.