GRB 201104A: A"Repetitive"Short Gamma-Ray Burst?

Gamma-ray bursts (GRBs) are divided into short gamma-ray bursts (SGRBs) and long gamma-ray bursts (LGRBs) based on the bimodal distribution of their durations. LGRBs and SGRBs are typically characterized by different statistical characteristics. Nevertheless, there are some samples that challenge such a framework, such as GRB 060614, a long-duration burst with short-burst characteristics. Furthermore, GRBs are generally considered to be an event with no periodic or repetitive behavior, since the progenitors usually undergo destructive events, such as massive explosions or binary compact star mergers. In this work, we investigated Fermi data for possible quasiperiodic oscillations and repetitive behaviors of GRBs using timing analysis methods and report a special event GRB 201104A, which is a long-duration burst with the characteristics of an SGRB, and it exhibits a"repetitive"behavior. We propose that such a situation may arise from lensed SGRBs and attempt to verify it by Bayesian inference. In addition, we extend the spectral analysis to Bayesian inference. In spite of the existence of at least two distinct time periods with a nearly identical spectrum, there is no strong evidence that they result from a lensing GRB. Taking the gravitational-lensing scenario out of consideration, a long burst would resemble a short burst in its repetitive behavior, which presents a challenge for the current classification scheme.


INTRODUCTION
The gamma-ray bursts (GRBs) is the brightest explosion in the universe and is characterized by high variability. GRBs can be divided into short gamma-ray bursts (SGRBs) and long gamma-ray bursts (LGRBs) based on the bimodal distribution of the observed durations in the BATSE era (Kouveliotou et al. 1993).
LGRBs originate in star-forming regions in galaxies and are observed in association with massive star-collapse supernovae (Woosley 1993;Fruchter et al. 2006). Short bursts are located in low star-forming regions of their host galaxies, and they are believed to originate from compact binaries (Eichler et al. 1989;Narayan et al. 1992;Gehrels et al. 2005;Fong et al. 2009;Leibler & Berger 2010;Fong & Berger 2013;Berger 2014;Tunni-cliffe et al. 2014;O'Connor et al. 2022). The duration of the burst (T 90 ) is often inconsistent across energy ranges and different instruments (Qin et al. 2012;Bromberg et al. 2013), so other characteristics of GRBs are important to consider. For example, the delay in the arrival of photons of different energies (Norris et al. 2000;Norris & Bonnell 2006), the E γ,iso and E p,z correlation (Amati et al. 2002), as well as the star formation rate of the host galaxy (Li et al. 2016). in the spectral analysis results of GRBs observed by Fermi, the low-energy spectral index of long bursts is lower than that of short bursts (Nava et al. 2011). However, there are still some special cases, such as long-duration bursts with short-burst characteristics (GRB 060614 (Gehrels et al. 2006;Yang et al. 2015;Jin et al. 2015) and GRB 211211A (Rastinejad et al. 2022;Gompertz et al. 2022;Yang et al. 2022)) and short-duration bursts with long burst characteristics (GRB 040924 (Fan et al. 2005) and GRB 200826A (Zhang et al. 2021;Ahumada et al. 2021;Rossi et al. 2022)).
All these issues are intertwined, so it is challenging to interpret all observed properties in a consistent picture. The confusion stems from our lack of understanding of the central engine and dissipation mechanisms. Studying the power density spectra of the prompt emission of GRBs can help to solve these puzzles (Beloborodov et al. 1998;Dichiara et al. 2013a). While most observations and theories indicate that GRBs do not repeat (Meszaros 2006;Kumar & Zhang 2015), some work has attempted to explore the possibility of quasi-periodic oscillation (QPO) signals in GRBs (Dichiara et al. 2013b;Zhang et al. 2016;Tarnopolski & Marchenko 2021). However, some bursts with repetitive behavior on the light curve are thought to be caused by gravitational lensing (Paynter et al. 2021;Wang et al. 2021;Veres et al. 2021;Lin et al. 2021). In this work, we investigated possible QPOs and repetitive behaviors in the prompt emission of GRBs observed by Fermi. We use the Bayesian test to find and confirm periodic or quasiperiodic oscillations in red noise, which builds upon the procedure outlined by Vaughan (2010). To identify possible repetitive events, we calculate the autocorrelation function proposed by Paynter et al. (2021). We have found a special example, GRB 201104A, and we have taken into consideration the possible repetition artifacts caused by the lens effect. Furthermore, we extend the spectral analysis to Bayesian inference of lens and nonlens models. We believe that this procedure should be taken into account in future research to certify gravitationally lensed GRBs.
The paper is organized as follows: In Section 2, we describe the search method we use to identify repeating events and QPOs in Fermi's GRBs. In Section 3, we present observations of GRB 201104A and a detailed analysis of its properties. In Section 4, we performed Bayesian inferences under the lensing hypothesis, considering both light curves and spectrum data. In Section 5, we summarize our results with some discussion.

COMPREHENSIVE SEARCH
In modern astrophysics, the analysis of time series is an essential tool. These two methods were used to search for possible QPO and repetitive behaviors in the current Fermi-GBM data (Meegan et al. 2009).

Autocorrelation function
Signal autocorrelation can be used to measure the time delay between two temporally overlapping signals. It may be an intrinsic property of the source, or it may be due to gravitational lensing. The standard autocor-relation function (ACF) is defined as follows: (1) To fit the ACF sequence, we apply the Savitzky-Golay filter F (δt). The values of the window length and the order of the polynomial are set to be 101 and 3, respectively (Paynter et al. 2021). The dispersion (σ) between the ACF and the fit F (k) is where N is the total number of bins. As usual, we identify the 3σ outliers as our candidates.

Power Density Spectra
To identify possible quasiperiods in the time series data, red noise must be modeled in order to assess the significance of any possible period. For the above purpose, we developed a procedure based on Vaughan (2005,2010,2013) and Covino et al. (2019). Also, we refer to Beloborodov et al. (1998);Guidorzi et al. (2012Guidorzi et al. ( , 2016; Dichiara et al. (2013a,b) for details on analyzing power density spectra in GRBs. Power density spectra (PDS) are derived by discrete Fourier transformation and normalized according to Leahy et al. (1983).
We consider two models to fit the PDS (Guidorzi et al. 2016), the first is a single power-law function plus white noise called a power-law model, here N is a normalization factor, f is the sampling frequency, and its lower limit is related to the length of the time series, which is 1/T (the time interval). The upper limit of f is the Nyquist frequency, which is 1/(2δ t ), where δ t is the time bin size of data. The value of white noise B is expected to be two, which is the expected value of a χ 2 2 distribution for pure Poissonian variance in the Leahy normalization. In some GRBs, there is a distinct broken power law. As such, we consider another model called BPL, There is one more parameter here, the break frequency  (Vaughan 2010). After we have the posterior distribution of the model parameters, we calculate the global significance of every frequency in the PDS according to T R = max j R j , where R = 2P/S, P is the simulated or observed PDS, and S is the best-fit PDS model. This method selects the maximum deviation from the continuum for each simulated PDS. The observed T R values are compared to the simulated distribution and significance is assessed directly. The corrections for the multiple trials performed were included in the analysis because the same procedure was applied to the simulated as well as to the real data.

OBSERVATION AND DATA ANALYSIS
We analyzed the GRBs in the Fermi GBM Burst Catalog (Gruber et al. 2014;Bhat et al. 2016;Von Kienlin et al. 2014, 2020 using the above method. In total, we examined the PDS of 248 short bursts and 920 long bursts (bounded by T 90 of 3 seconds) with peak counts greater than 50 in 64 ms time bin. Among these samples, we did not find a quasi-periodic signal with a significance exceeding 3 σ, but among the candidates of the ACF check, we found an interesting example,which is GRB 201104A. The basic observation information and data reduction are as follows.
The Fermi-GBM team reports the detection of a possible long burst GRB 201104A (trigger 626140861.749996 / 201104001) (Fermi GBM Team 2020). At the same time, Fermi-LAT (Atwood et al. 2009) also observed high-energy high photons from this source with high significance (Ohno et al. 2020). In addition, we try to search the observation data of other telescopes through GCN 1 , but there is no X-ray and optical data available with other telescopes.
We present a further analysis of GRB 201104A observed by Fermi instruments in this work. Fermi-GBM (Meegan et al. 2009) has 12 sodium iodide (NaI) detectors and 2 bismuth germanate (BGO) detectors. According to the pointing direction and count rate of the detectors, we selected a NaI detector (n8) and a BGO detector (b1) respectively. The Fermi-GBM data is processed with GBM Data Tools (Goldstein et al. 2021), which makes it incredibly easy for users to customize. We performed a standard unbinned maximum likelihood analysis for GRB 2 using Fermitools 3 (version 2.0.8), and we determined the probability of each photon originating from this source. In Figure 1 (a), we present the GBM and LAT light curves for several energy bands. It appears to be a long burst of three episodes and with a weak extending component. ACF and PDS results are also shown in (b) and (c) in Figure 1. A double peak in the ACF means that the light curve is similar after two time shifts, while there is also a peak in the PSD at a frequency of 0.3 Hz. While the significance level of this quasi-periodicity is not high enough, we also analyze in detail each possible episode in the following section.

Spectral Analysis
In order to confirm the duration of this burst, we recalculated the T 90 (Koshut et al. 1996) of the GBM n8 detector in the energy range of 50-300 keV. We then used the Bayesian block technique (Scargle et al. 2013) to determine the time interval of this burst (see the left of Figure 2). We divide the main burst [T 0 -0.1 , T 0 + 8.3 s] into three episodes, which are Episode a [T 0 -0.1 , T 0 + 2.7 s], Episode b [T 0 + 2.7, T 0 + 5.5 s], and Episode c [T 0 + 5.5, T 0 + 8.3 s]. We perform both time-integrated and time-resolved spectral analyses of GRB 201104A, and the specific time interval is shown in Table 1. For each time interval, we extract the corresponding source spectra, background spectra and instrumental response files following the procedure described Section 4.3 and Goldstein et al. (2021). In general, the energy spectrum of GRB can be fitted by an empirical smoothly joined broken power-law function (the so-called "Band" function; Band et al. 1993). The Band function takes the form of Where A is the normalization constant, E is the energy in unit of keV, α is the low-energy photon spectral index, β is the high-energy photon spectral index, and E 0 is the break energy in the spectrum. The peak energy in the νF ν spectrum is called E p , which is equal to E 0 ×(α+2). In addition, when the count rate of high-energy photons is relatively low, the high-energy spectral index β often cannot be constrained. In this case, we consider using a cutoff power-law (CPL) function, where α is the power law photon spectral index, E c is the break energy in the spectrum, and the peak energy E p is equal to E c × (2 + α).
In our joint spectral fitting analysis of Fermi-GBM and Fermi-LAT, we use the pgstat statistic 4 in Bayesian inference. The process used is the same as that described in Section 2.2. The difference is in the model and likeli-hood function, as well as the additional process of folding the model and the instrumental response file. The fitting results of the spectrum are presented in Table  1, and the evolution of the time-resolved spectrum is illustrated in Figure 2 (b). Evidently, the spectral evolution does not quite satisfy the hard-to-soft evolution pattern or intensity-tracking pattern (Norris et al. 1986;Golenetskii et al. 1983;Lu et al. 2012). Furthermore, the posterior parameter of the spectra analysis of Episode a and Episode b are very similar, which will be further analyzed in Section 4.

Empirical Correlations and T 90 Related Distribution
Through the spectral analysis result in Section 3.1, we attempt to classify GRB 201104A using the Amati correlation (Amati et al. 2002;Minaev & Pozanenko 2020).Due to the lack of a measured redshift, we calculated E γ,iso and E p,z in different redshifts (from 0.1 to 5) with the cosmological parameters of H 0 = 69.6 kms −1 Mpc −1 , Ω m = 0.29, and Ω Λ = 0.71. The star signs of the different colors in Figure 3 (a) display the results obtained with different redshifts for each episode, and clearly they all depart from the range of LGRBs (type II).
Additionally, Minaev & Pozanenko (2020) proposed a new classification scheme combining the correlation of E γ,iso and E p,z and the bimodal distribution of T 90 . The parameter EH is proposed to characterize the Amati correlation, The T 90,z -EH trajectories calculated in different redshifts (from 0.01 to 5) for each episode of GRB 201104A are shown in Figure 3 (b). In the calculation of the above two correlations, we made assumptions about the redshift. For the Amati correlation, within our hypothetical redshift range, each episode is an outlier beyond the 2σ of the LGRBs classification. In the T 90,z -EH correlation, each epsiode is also within the region of the SGRBs. But for the main burst, Its trajectory will cross the boundary line between SGRBs (type I) and LGRBs (type II). When the redshift is 0.26 < z < 3.42, it is located in the region of the LGRBs. For the T 90 episode, a higher redshift is required to fall into the classification of SGRBs. However, in some current work (Paterson et al. 2020;O'Connor et al. 2022;Hjorth et al. 2012;Vergani et al. 2015;Palmerio et al. 2019), the redshift of the SGRBs category falls between 0.5 and 2. Therefore, we propose a possible situation to make the short burst duration longer in Section 4.
We also examine other T 90 -related distributions through the bursts list given by Fermi GBM Burst catalog (Gruber et al. 2014;Bhat et al. 2016;Von Kienlin et al. 2014, 2020. We collected the E p and fluence of each burst and calculated the hardness ratio (HR), which is the ratio of the observed counts in 50-300 keV compared to the counts in the 10-50 keV band (Goldstein et al. 2017). In Figure 3 (c), (d) and (e), we plotted the HR, E p and fluence of each episode for GRB 201104A and other catalog's bursts together, and fitting the distribution by using a two-component Gaussian mixture model with scikit-learn (Pedregosa et al. 2011). When the extended emission is not considered, each episode will gradually approach SGRB (type I) in classification.

Spectral Lag
In most GRBs, there is a lag between different energy bands, which is called the spectral lag. In general, LGRBs exhibit a relatively significant spectral delay (Norris et al. 2000;Gehrels et al. 2006), whereas short hard bursts display a negligible spectral lag (Norris & Bonnell 2006). Besides, a fraction of short GRBs even show negative lags (Yi et al. 2006). A cross-correlation function (CCF) can be used to quantify such an effect since the pulse peaks at different energy bands are delayed. Unlike ACF, CCF is used to compare different time series. The method is widely used to calculate spectral lag (Band 1997;Ukwatta et al. 2010). We calculated the CCF function for GRB 201104A in the energy bands between 100-150 keV and 200-250 keV from T 0 -0.1 to T 0 +8.3 s (the main burst episode), and calculated the peak value of CCF after polynomial fitting. By using Monte Carlo simulations, we can estimate the uncertainty of lags (Ukwatta et al. 2010). The spectral lag of GRB 201104A is (τ = -28.0 ± 56.1) ms and is compared with the delays of other long and short GRBs (Bernardini et al. 2015), as shown in Figure 4. In the observer frame, the mean values of spectral lags for long and short GRBs are (τ L = 102.2 ± 38.1) ms and (τ S = -0.73 ± 7.14) ms (Bernardini et al. 2015), respectively. Although the result of the spectral lag of this burst is closer to the mean value of SGRBs, it cannot be classified accurately.

LENSING HYPOTHESIS
To explain the similar spectra of Episode a and Episode b and their long-duration but short-burst characteristics, we propose that Episode b is the result of lensing in Episode a and Episode c is a relatively soft and weak extended radiation, so no repeats of it have been observed.
In a gravitational lensing system, photons that travel longer distances arrive first, because a shorter path means that the light passes through deeper gravitational potential well of the lens, where the time dilation is stronger. The source flux is lower for photons arriving later compared to those which arrive earlier. There will therefore be at least one early pulse followed by a weaker pulse for a lensed GRB. The time delay between these two pulses is determined by the mass of the gravitational lens. For lensing of a point mass, we have (Krauss & Small 1991;Narayan & Wallington 1992;Mao 1992) ( where ∆t is the time delay, r is the ratio of the fluxes of the two pulses, and (1+z l )M l is the redshifted lens mass. With the measured ∆t and r, it is straightforward to calculate the redshifted mass (1+z l )M l . Using Bayesian inference methods of light curve and energy spectrum data, we will estimate the parameters and compare lens and non-lens models.

Light-curve Inference
Paynter et al. (2021) developed a Python package called PyGRB to create light-curves from either prebinned data or time-tagged photon-event data. We extend the analysis to the Fermi data as well (Wang et al. 2021). Here the same method as in Section 2.2 is used to obtain the posterior distributions of the parameters. Bayesian evidence (Z) is calculated for model selection and can be expressed as where θ is the model parameters, and π(θ) is the prior probability. For TTE data from various instruments, the photon counting obeys a Poisson process and the likelihood ln L for Bayesian inference takes the form of where N i stands for observed photon count in each time bin, and the model predicted photon count consists of the background count δt i B and the signal count δt i S(t i |θ). Note that the differences of Z among models are important for our purpose. Usually the light curve of a GRB is a pulse shape of fast-rising exponential decay (FRED), where ∆ is the start time of pulse, A is the amplitude factor, τ is the duration parameter of pulse, and ξ is the asymmetry parameter used to adjust the skewness of the pulse. Through different FRED functions, we define different light curve models S(t i |θ) to describe whether the pulses are lensed images or not, the lensing and null scenarios as S non-lens (t|θ non-lens ) =S(t|∆ 1 , A 1 , τ 1 , ξ 1 ) For lens model, r is the flux ratio between two pulses (see Equation (8)) and B is a constant background parameter. The ratio of the Z for two different models is called as the Bayes factor (BF) and the logarithm of the Bayes factor reads ln BF lens non-lens = ln(Z lens ) − ln(Z non-lens ).
As a statistically rigorous measure for model selection, if ln (BF) > 8 we have the "strong evidence" in favor of one hypothesis over the other (Thrane & Talbot 2019). This method was used to analyze the time series of Episode a and Episode b. We masked the light curve after Episode b with Poisson background in order to exclude the influence of other time periods. In Table  2, the results of Bayesian inference based on the light curves of the NaI and BGO detectors are presented.

Pearson Correlation Coefficient
The pulse shape we observe in reality does not correspond to a simple FRED function, and if multiple FRED functions are used to construct the model, Bayesian inference will become significantly more difficult. Calculating the Pearson correlation coefficient is a simple and effective method of time series analysis. Pearson correlation coefficients were calculated for Episode a and Episode b for two different energy bands (see Figure 5). In the NaI detector energy band (50-300 keV), the results are: r = 0.61, p = 1.24 × 10 −5 , and in the BGO detector energy band (300-40000 keV), the results are: r = 0.43, p = 3.80 × 10 −3 . The results suggest that there is a general correlation between Episodes a and b.

Spectral Inference
Under the lensed GRB hypothesis, in addition to requiring the same shape of the light curves, the similarity of the spectrum must also be considered. The cumulative hardness comparison in different energy bands is a simple but statistically powerful methodology (Mukherjee & Nemiroff 2021). And such a method has been used in some works (Wang et al. 2021;Veres et al. 2021;Lin et al. 2021) as one of the indicators to confirm the lensed GRB.
We propose a procedure that considers spectral fitting with Bayesian inference for lens and non-lens models in order to achieve this goal. Typically, the detector response files of GBM contain one or more detector response matrices (DRMs) encoding the energy dispersion and calibration of incoming photons at various energies to recorded energy channels (Goldstein et al. 2021). The matrix also encodes the effective area of the detector as a function of energy relative to the source to which the detector is pointing. Due to the strong angular dependence of the response (and the degree of angular dependence varies with energy), the effective area can fluctuate significantly. Therefore, we should select the DRM that corresponds to the period we are interested in. By interpolating the time series DRMs, we can obtain the DRM of the center time of the interaction we are interested in, and generate the corresponding response file. It is relatively simple to get the response file of LAT, which is generated by using gtrepden in Fermitools. Since GRBs are relatively short in duration, so that the accumulation of LAT background is negligible. To obtain the GBM background file, we use polynomial fitting for the two periods before and after the burst, and then interpolate to obtain the background of the selected time. However, although background and instrument responses do not differ significantly in the same GRB event, it is necessary to consider these variations when searching for lensing effects in different GRB events.
The likelihood function used in the spectral inference is the pgstat mentioned in Section 3.1. And we use a Band function (Equation 5) and a ratio parameter r to construct the lens model, The non-lens model is composed of two independent Band functions, N non-lens (E|θ non-lens ) =N Band1 (E|α 1 , β 1 , E c1 , A 1 ) + N Band2 (E|α 2 , β 2 , E c2 , A 2 ).
It should be noted that N band1 and N band2 will fold the response files of Episode a and Episode b respectively. We use the above model to fit the spectra of the two episodes. The method of model comparison is consistent with Section 4.1. As shown in Table 2, the results of the spectral inference (ln (BF) = 8.21) are more inclined to suggest that the spectra of these two time periods differ only in the normalization constant (see Figure 6). In fact, the spectral evolution of GRBs is very significant and common, and E p -Flux usually has two types of evolution patterns, hard-to-soft evolution pattern and intensity tracking pattern (Norris et al. 1986;Golenetskii et al. 1983;Lu et al. 2012), respectively. So for GRB 220411A, the nearly identical spectrum of the first two episodes is peculiar.

SUMMARY AND DISCUSSION
In this work, we investigated Fermi's GRBs for possible QPOs and repetitive behavior events, as well as performed a detailed analysis of burst 201104A. Our findings are the following: • Following the current analysis method in this work, there is no significant QPO signal greater than 3 σ in the light curve of 248 short bursts and 920 long bursts selected. However, some GRBs exhibit repetitive behavior, such as GRB 201104A.
• GRB 201104A has similar temporal evolutions for both Episodes a and b, as well as the posterior parameters of the spectral fitting result. According to the Amati correlation diagram, each episode is closer to the group of SGRBs. Classifying it as a short burst is also favored in the T 90 -related distributions when the extended emission component is not taken into account. Consequently, GRB 201104A is a long-duration burst with characteristics of SGRBs.
• The Bayesian inferences of the light curve do not fully support the lens model. Nevertheless, the spectral inference result supports the lens model, at least showing that the spectra of these two episodes are very consistent.
A long-duration burst with the characteristics of a short burst can be explained very naturally by the gravitational-lensing scenario, although this is a very coincidental circumstance. As well, it is not impossible that this event may be just a very special occurrence that comprises at least two intrinsically similar episodes. In other words, there is a repetition mechanism in GRB, like the central engine with memory that Zhang et al. (2016) once proposed. This type of burst has many features characteristic of a SGRB, but unless the Li-Paczynski macronova (also known as the kilonova) can be detected (Li & Paczyński 1998), this will undoubtedly be the strongest evidence for a SGRB of compact-binary origin, such as GRB 060614 (Gehrels et al. 2006;Yang et al. 2015;Jin et al. 2015) and GRB 211211A (Rastinejad et al. 2022;Gompertz et al. 2022;Yang et al. 2022).
Furthermore, our proposed procedure is essentially consistent with the method proposed by Mukherjee & Nemiroff (2021) for comparing the cumulative hardness in different energy bands. In this procedure, we use the Bayes factor to determine the consistency of the spectrum. Over time series analysis, the advantage is that low count rate instruments can be considered in Bayesian inference, such as the spectra of high-energy photons detected by Fermi-LAT. For certifying different GRB events as gravitational lensing events, this procedure could take into account the effects of instrumental responses as well as background. There is currently a problem with this procedure, which is that the episodes used for spectral inference must be selected in advance. Therefore, this procedure will only give a posterior dis-tribution of the ratio r without a time delay caused by gravitational lensing. We can set the time delay and time window as free parameters to solve this problem. Since each model calculation must take into account the change in the spectral fitting file, the calculation time will greatly increase. Our future work will optimize this step in order to search for gravitationally lensed GRBs. Currently, the most complete procedure is to use both time series and spectrum data to perform Bayesian inference for lens and non-lens models.

ACKNOWLEDGMENTS
We appreciate the anonymous referee for their helpful suggestions. We thank S. Covino and Yi-Zhong Fan for their important help in this work. We appreciate Zi-Min Zhou for his helpful suggestions. We acknowledge the use of the Fermi archive's public data. This work is supported by NSFC under grant No. 11921003.