Signatures of the Self-organized Criticality Phenomenon in Precursors of Gamma-ray bursts

In a fraction of gamma-ray bursts (GRBs), some weaker pulses prior to main bursts have been detected, which are usually called as precursors. Single and multiple precursors have been found in both short GRBs and long GRBs (e.g. Troja et al., 2010; Hu et al., 2014; Li & Mao, 2022). Some studies reported that there is no obvious difference between precursors and main bursts by the temporal and spectral analysis, indicating that they might originate from the same mechanism and share the same origin (e.g. Burlon et al., 2008, 2009; Lan et al., 2018; Li et al., 2021, 2022). Some theoretical models have been proposed to explain the origin of precursors, involving photospheric precursors (e.g., Murakami et al., 1991; Lyutikov & Usov, 2000; Daigne & Mochkovitch, 2002), shock breakout precursors (e.g., Ramirez-Ruiz et al., 2002; Wang & Mészáros, 2007) and so on. It is expected that precursors provide important clues as to the nature of the central engines of GRBs, and more importantly, unveil the physical mechanism of themselves, which has so far remained mysterious.

For a self-organized criticality (SOC) system, subsystem with a nonlinear energy dissipation process in an external field assembles itself into a critical state. Owing to a small perturbation, the critical state is triggered and will generate an avalanche-like chain reaction of any size (Bak et al., 1987). Typical characteristics of an SOC system are that the event sizes present scale-free power-law-like distributions (Aschwanden, 2011; Wang & Dai, 2013; Yi et al., 2016; Zhang et al., 2022). A universal statistical analytical SOC model developed by Aschwanden (2012) provided an theoretical framework to quantitatively connect the concept of fractal dimensions to the power-law distributions. In addition, another remarkable hallmark of the SOC system is the scale-invariance structure of the avalanche size differences, which can be described by the probability density functions (PDFs) with a invariant $q$-Gaussian form (Caruso et al., 2007; Wei et al., 2021; Wang et al., 2023).

Aschwanden (2022) reported that the dynamical nonlinear behaviors of almost all of the discipline in astrophysical systems are governed by SOC process. As the most luminous explosion events in the universe, the physical processes and radiation mechanisms of GRBs are still debated. The sizes frequency distributions of GRB X-ray flares, including the energies, durations, peak fluxes, rise times, decay times and waiting times, etc., have been reported to be attributed to a magnetic reconnection process on the basis of the fractal-diffusive avalanche (Wang & Dai, 2013; Yi et al., 2016; Wei, 2023). Besides X-ray flares, similar size distributions of GRB optical flares are found by Yi et al. (2017), indicating the similar SOC framework to X-ray flares. It is generally believed that X-ray flares may be internal emission through the collisions of relativistic shells ejected due to the delayed central engine activities (e.g. Burrows et al., 2005; Zhang, 2006; Mészáros, 2006). Thus, it is necessary to investigate whether the GRB prompt phases can also be attributed to an SOC process. It is worth noting that the SOC behavior in GRB prompt phase have been studied by Lyu et al. (2020) and the results tentatively suggested that SOC phenomena may exist. However, only the GRBs with at least three pulses have been investigated yet. In practice, a fraction of short GRBs are overwhelmingly single-peaked and double-peaked bursts (Hakkila et al., 2018; Li et al., 2020, 2021), hence we have examined in detail the properties of the single-peaked and double-peaked BATSE and Swift short GRBs and confirmed the SOC behaviors in short GRB prompt phases in our recent work (Li et al., 2023, hereafter Paper I).

As an important component of GRB prompt phase, precursor provides important information on the nature of central engines. However, the indisputable fact is that no work has been done yet to disclose whether the SOC features also exists in GRB precursors. Therefore, we perform a analysis to search for the evidences of the SOC behaviors in precursors in this letter. We focus on the size distributions and further check the scale-invariance structure of the size differences. In addition, we also investigate the two kinds of distributions of main bursts to examine the possible connection between them. This paper is organized as follows. Data analysis and methods are given in Section 2. Our main results are presented in Section 3, in which we will give special attention to a direct comparison with the results of Paper I. Finally, the discussions and conclusions are presented in Sections 4 and 5.

Our precursor sample is from Li & Mao (2022). They limited the bursts to events with durations longer than 2 s from 1006 GRBs in the third Swift/BAT catalog. The signature of the precursor sample should be clearly shown in all Swift energy-channel 64 ms binning lightcurves, including 15 - 25, 25 - 50, 50 - 100, 100 - 350, and 15 - 350 keV. Finally, they identified 52 long GRBs with the precursor activity. In total, they obtained 227 pulses, including 127 main pulses and 100 precursor pulses. Then, they use the fitting function provided by Norris et al. (2005) to fit all the precursor and the main burst pulses. In addition, peak-time intervals and quiescent times between different episodes were extracted. The detailed sample information and data processing refers to Li & Mao (2022). In this letter, we select the characteristic parameters of the precursor pulses to analyse, including peak flux ($P$), duration ($D$), rise time ($t_{\rm r}$), and decay time ($t_{\rm d}$). Besides these parameters, we also investigate 52 peak-time intervals ($\Delta t_{\rm p1}$) and 52 quiescent times ($\Delta t_{\rm p2}$) between precursors and main bursts. Note that considering the limited number of GRBs with more than two precursors from Li & Mao (2022), we thus only select the peak-time intervals and quiescent times of GRBs with one precursor.

As applied in our Paper I, an thresholded size distribution function motivated by physical arguments (e.g. instability threshold, detection threshold, background noise) (Aschwanden, 2015) will be utilized again, which can be written as

where $N$ is the number of events, $\alpha_{\rm d}$ is the power-law index, $x_{\rm 1}$ represents the minimum cutoff values of scale-free range, $x_{\rm 2}$ is the upper bound by the largest event, which represents a sharp cutoff in the differential size distribution, and $x_{\rm{0d}}$ is a threshold value of differential size distribution, which is common in astrophysical data sets and represents the deviations from ideal power law size distributions (Newman, 2005) caused by some natural effects. $\sigma_{d}=\sqrt{N_{bin,i}}/\triangle x_{i}$ represents the uncertainty of the differential size distribution, where $N_{bin,i}$ and $\triangle x_{i}$ are the number of events of the $i$-bin and the bin size, respectively. Note that such the thresholded power law size distribution is also called a “Generalized Pareto distribution”, the “Generalized Pareto Type II distribution” and the “Lomax distribution” (e.g. Hosking & Wallis, 1987; Johnson, 1994; Lomax, 1954). For small data sets and near the upper cutoff, it is statistically advantageous to use the cumulative size distribution of the Equation 1 (Aschwanden, 2015), which includes all events accumulated above some size $x$ and can be obtained as

where $N_{ev}$ refers to the total number of events, $\alpha_{\rm c}\neq 1$ is the power-law index, and $x_{\rm{0c}}$ is a threshold value of cumulative size distribution. The uncertainty of the cumulative distribution in a given bin $i$ is estimated with $\sigma_{c}=\sqrt{N_{i}}$, where $N_{\rm i}$ is the number of events of the bin. Due to the small sample, a rank-order plot which is essentially an optimum adjustment to minimum statistics is used to get the size distributions. The detailed method refers to the related descriptions in Aschwanden (2011, 2019).

The Markov chain Monte Carlo (MCMC) technique is first utilized to obtain the best-fit results with the package $pymc$. Then, the standard reduced chi-square ($\chi{{}_{\rm\nu}^{2}}$) criterion which is more suitable for modeling astrophysical data with some natural effects is used to test the goodness of fits (Aschwanden, 2019), although it may give biased and noisy slope estimates (Goldstein et al., 2004; Bauke, 2007). The $\chi{{}_{\rm\nu}^{2}}$ criterions can be written as

for the differential distribution function, and

for the cumulative distribution function (Aschwanden, 2015), where $n_{\rm x}$ is the number of logarithmic bins, $n_{\rm p}$ is the number of the free parameters, $N_{\rm fit,d}(x_{i})$ and $N_{\rm fit,c}(x_{i})$ are the theoretical values, and $N_{\rm obs,d}(x_{i})$ and $N_{\rm obs,c}(x_{i})$ are the corresponding observed values, respectively. Note that the points below the threshold $x_{0}$ are just noise and do not contribute to the accuracy of the best-fit power-law index, thus they are ignored when the reduced chi-square is calculated.

The avalanche size difference is defined as $W_{n}=Z_{i+n}-Z_{i}$, where $n$ is the temporal interval scale and $Z_{i}$ is the scale size of the $i{th}$ burst. Here the dimensionless difference is rescaled to $w_{n}=W_{n}/\sigma_{W_{n}}$, where $\sigma_{W_{n}}$ is the standard deviation of $W_{n}$. The Tsallis q-Gaussian function given by Tsallis et al. (1998) is used to fit the PDFs of the size differences, which can be written as

where $q$, $a$, and $b$ are free parameters and correspond to the sharpness of the peak, the normalization factor, and the width of the peak, respectively. Note that $q$ describes the departure from the Gaussian function. In other words, when $q\rightarrow 1$, the $q$-Gaussian form reduces to the Gaussian function with mean $\mu=0$ and standard deviation $\sigma=1/\sqrt[]{2b}$. Here, due to the limited number of data points, we use the cumulative distribution function (CDF) of Tsallis $q$-Gaussian function to fit the size difference distributions (Sang & Lin, 2022; Wang et al., 2023), which can be written as

The best-fit results are obtained by minimizing the $\chi^{2}$ statistics.

According to the prediction proposed by Aschwanden (2012, 2014) for a fractal-diffusive SOC model, the indexes $\alpha_{P}$ and $\alpha_{T}$ of the power-law size distributions for the peak flux and duration can be derived from $\alpha_{P}=2-1/S$ and $\alpha_{T}=(1+S)\beta/2$, where $S=1,2,3$ are the spatial dimensions of the SOC model and $\beta=1$ is the diffusive spreading exponent for the classical diffusion. Obviously, the spatial dimension $S=3$ can be obtained to explain the SOC behavior of precursors by taking the mean values of the power-law indexes. Likewise, the same spatial dimension $S=3$ can be determined by the quite same mean values of the power-law indexes of the main bursts, indicating that precursors and main bursts might share the similar physical processes. Moreover, we find that the results of main bursts are consistent with those reported by our Paper I (Li et al., 2023) and Lyu et al. (2020) in despite of different GRB samples, involving short or long or mixed GRBs samples, suggesting that the similar processes are responsible for producing the main bursts in two classes of GRBs. In particular, under the assumption that there is no correlation between the sizes of two events, an exact relation between the index $\alpha$ of the avalanche size distribution and the $q$ value of the appropriate $q$-Gaussian function has been obtained as (Caruso et al., 2007; Celikoglu et al., 2010; Wei et al., 2021)

Thus, the theoretical $q$ values can be derived as $q_{P}=2.2$ and $q_{T}=2$ by taking the theoretical values of $\alpha$ when $S=3$. In fact, we find that the $q$ values of precursors are found to be slightly larger than those of main bursts. However, the $q$ values exhibited by both the precursors and the main bursts are smaller than the theoretical ones.

At present, the fact remains that whether precursors are inherently different from main bursts because there is no unambiguous definition of “precursor”. One is that precursors may indeed be generated from the different physical processes. The other is that precursors possess the similar origins to those of the main bursts and the absence of the dimmer precursors might be due to the limited sensitivity or energy coverage of the current GRB detectors (e.g. Li et al., 2021). In the fireball model, the central engines of GRBs have long-lasting activities and arise the prompt emissions from early internal shocks and X-ray flares from late internal shocks through collisions of relativistic shells ejected. Precursors could be photospheric emission when it is bright enough to be detected as a separate component in time from the main bursts (e.g. Lyutikov & Usov, 2000). In this scenario, the precursors are favorable to the non-thermal emission, which are not well consistent with the previous spectral analysis results (Burlon et al., 2008; Hu et al., 2014; Li et al., 2022). The mean value of quiescent time in the precursor sample is 14.6 s and the longest one is $\sim$ 526 s (Li & Mao, 2022), which can be explained by multistage collapse model (Lipunova et al., 2009) and collapsar fallback model (Wang & Mészáros, 2007). In addition, the GRB jet may be magnetically dominated and the magnetic energy can be dissipated by magnetic reconnection based on the quiescent period timescale (Lyutikov & Blandford, 2003). Dănilă et al. (2015) reported that GRBs can be produced in a cellular automaton model of a magnetic flux tube with background flow and explained in an SOC framework via magnetic reconnection processes in the magnetized fireball. Moreover, they found that the energetic initial impulses in the magnetohydrodynamic flow can cause one-dimensional signatures when a two-dimensional SOC regime has been established. Owing to the radial component of magnetic fields decaying with the radius, the early prompt phases are in a three-dimensional form while the later X-ray flares are closer to one dimension (Wang & Dai, 2013; Lyu et al., 2020; Li et al., 2023). Thus, the same Euclidean dimensions of SOC systems between precursor and main burst near to the central engine is reasonable.

Recently, a peculiar long GRB 211211A is found to be associated with a kilonova, implying that the event may originate from a merger of compact star binaries (e.g. Yang et al., 2022; Troja et al., 2022). About 1 s prior to the main burst, a precursor flare lasting for $\sim$ 0.2 s presents a feature of quasi-periodic oscillation at frequency $\approx$ 22 Hz, suggesting that the flare results from a resonant shattering of a magnetar before coalescence (Xiao et al., 2022; Suvorov et al., 2022). On the other hand, the detection of fast radio burst (FRB) 200428 from the Galactic soft gamma-ray repeater (SGR) J1935+2154 firstly confirms the magnetar origin of some FRBs (CHIME/FRB Collaboration et al., 2020; Bochenek et al., 2020). In addition, sub-second periodicity in FRB 20191221A further strengthens the evidence of magnetar origin (Chime/Frb Collaboration et al., 2022). It is worth noting that the power-law size distribution features and scale-invariance structures in FRBs have been reported, indicating that both FRBs and GRBs can be explained by the same physical framework of fractal-diffusive SOC systems with the spatial dimension $S=3$ (Wang & Yu, 2017; Wei et al., 2021; Zhang et al., 2021; Wang et al., 2023). Thus, it would be very intriguing to contain whether the central engines of at least some GRBs are magnetars by the SOC behaviors. Even the physical correlation between GRBs and FRBs can be further studied in the future.

In this letter, we concentrated on investigating the SOC signatures of precursors in long GRBs. We found that similar scale-free power-law distribution properties and scale-invariance structure of avalanche size differences all exist in precursors, which can be well explained with the SOC framework driven by a magnetically dominated process. In addition, we further inspected the SOC behaviors of main bursts and confirmed that precursors and main bursts have similar statistical properties, indicating that both of them follow likely the same stochastic process and physical mechanism. The similar power-law indexes of size distributions and $q$ values of $q$-Gaussian distributions can be attributed to the same SOC systems with spatial dimension $S=3$.

We would like to thank Huang Yong-Feng, Zhang Zhi-Bin and Yi Shuang-Xi for beneficial discussions and helpful suggestions. We greatly acknowledge the anonymous referees for the valuable comments. This work is supported by Shandong Provincial Natural Science Foundation (Grant Nos. ZR2023MA049 and ZR2021MA021).