A Search for Millilensing Gamma-Ray Bursts in the Observations of Fermi GBM

Millilensing of Gamma-Ray Bursts (GRBs) is expected to manifest as multiple emission episodes in a single triggered GRB with similar light-curve patterns and similar spectrum properties. Identifying such lensed GRBs could help improve constraints on the abundance of compact dark matter. Here we present a systemic search for millilensing among 3000 GRBs observed by the \textit{Fermi} GBM up to 2021 April. Eventually we find 4 interesting candidates by performing auto-correlation test, hardness test, and time-integrated/resolved spectrum test to the whole sample. GRB 081126A and GRB 090717A are ranked as the first class candidate based on their excellent performance both in temporal and spectrum analysis. GRB 081122A and GRB 110517B are ranked as the second class candidates (suspected candidates), mainly because their two emission episodes show clear deviations in part of the time-resolved spectrum or in the time-integrated spectrum. Considering a point mass model for the gravitational lens, our results suggest that the density parameter of lens objects with mass $M_{\rm L}\sim10^{6} M_{\odot}$ is larger than $1.5\times10^{-3}$.


INTRODUCTION
In general relativity, light from very distant sources would be deflected by intervening masses, which is called gravitational lensing effect. Depending on the positions of the source, lens and observer, and the mass and shape of the lens, gravitational lensing can manifest itself by making one source produce multiple images, or altering the source image through shearing and convergence (Schneider et al. 1992). Due to the rapid progress in time domain surveys, many searches have been carried out for gravitational lensing of explosive transients, because such systems have been widely proposed as promising cosmological and astrophysical probes (Oguri 2019, for a review). For instance, transients strongly lensed by intervening galaxies could be used to improve constraints on cosmological parameters such as the Hubble constant H 0 Li et al. 2018), difference of time delays between multiple images among different particles, energy, or messengers, could be used for testing fundamental physics from the propagation speed (Biesiada & Piórkowska 2009;Fan et al. 2017;Collett & Bacon 2017), statistical search results of lensed transients with different time scales could be applied to derive constraints on the abundance of compact dark matter in different mass ranges (Muñoz et al. 2016;Liao et al. 2020;Zhou et al. 2022).
As one of the most violent explosions in the Universe, Gamma-Ray Bursts (GRBs) are bright enough to be detected in high-redshift range up to at least z ∼ 10 (Tanvir et al. 2009;Salvaterra et al. 2009), so they have long been proposed as the most promising transients for searching for gravitational lensing effect (Paczynski 1987;Mao 1992;Li & Li 2014). Thanks to the successful operation of several dedicated detectors, e.g., the Burst And Transient Source Experiment (BATSE) on Compton Gamma Ray Observatory (CGRO) (Meegan et al. 1992), the Burst Alert Telescope (BAT) on the Neil Gehrels Swift Observatory (Gehrels et al. 2004; Barthelmy et al. 2005) and the Gamma-Ray Burst Mon-itor (GBM) on the Fermi Observatory (Meegan et al. 2009), ∼ 10 4 GRBs have been detected. Two kinds of searches for lensed GRBs have been widely carried out: • searching for independently triggered GRB pairs with similar light curve and spectrum, but with different flux and a small positional offset (Nemiroff et al. 1994;Veres et al. 2009;Davidson et al. 2011;Li & Li 2014;Hurley et al. 2019;Ahlgren & Larsson 2020) 1 . This is usually called macrolensing event, for which the lens candidates are galaxies or clusters of galaxies (with mass > 10 7 M ).
• searching for two emission episodes in a single triggered GRB with similar light-curve patterns and similar spectrum properties. This is usually called millilensing event, for which the lens candidates could be an intermediate-mass black hole (Paynter et al. 2021), compact dark matter Nemiroff et al. 1993), star cluster or population III star (Hirose et al. 2006), with mass 10 4 ∼ 10 7 M .
Up to now, all searches for macrolensing events have yielded null results. On the other hand, however, several candidates of millilensing events have been proposed. Paynter et al. (2021) claims the first convincing evidence for millilensing events in the light curve of BATSE GRB 950830 and thus claims the existence of intermediatemass black holes 2 . Later on, an increase of millilensing GRB candidates from the data of Fermi/GBM have been claimed (Kalantari et al. 2021;Wang et al. 2021;Yang et al. 2021;Veres et al. 2021). For instance, Wang et al. (2021) and Yang et al. (2021) made independent analysis for GRB 200716C and argued it showing millilensing signatures. Veres et al. (2021) made an exhaustive temporal and spectral analysis to claim that GRB 210812A shows strong evidence in favor of the millilensing effects.
Inspired by the claims of these individual cases, it is essential to operate systematically search for millilensing events in the entire Fermi sample, because 1) after a decade of operation, thousands of GRBs were detected by the Fermi satellite, new candidates may be found; 2) Fermi/GBM is almost full-time monitoring the whole sky, the millilensing event rate for Fermi GRBs could provide powerful probes of the small-scale structure of the Universe.
1 Recently, Chen et al. (2021) suggests that the multi-band afterglow data should also be considered in the search for lensed GRB events. 2 Note that other searches based on the data of BATSE find none well-proved candidates (Nemiroff et al. 1993(Nemiroff et al. , 1994(Nemiroff et al. , 2001Ougolnikov 2003;Hirose et al. 2006;).
Most recently, Kalantari et al. (2021) applied autocorrelation method to the Fermi sample and find one more millilensing candidate (GRB 090717A). But Mukherjee & Nemiroff (2021a) later argued that the light curves of two pulses in GRB 090717 differ at about 5σ confidence level, therefore GRB 090717 does not present a compelling example of gravitational lensing. Nevertheless, in the analysis of Kalantari et al. (2021), spectral properties for each source have not been taken into account (not even the hardness test). In existing case studies, it can be seen that spectral analysis plays a vital role in identifying the authenticity of candidates. For instance, Mukherjee & Nemiroff (2021b) found cumulative hardness discrepancies between the two pulses in GRB 950830, thus argued that the case for GRB 950830 involving a gravitational lens may well be considered intriguing -but should not be considered proven.
In the study of GRB 210812A, Veres et al. (2021) shows that time-resolved spectrum is a more powerful tool than hardness to help justify the lensing effect. Taking advantage of its two main instruments (GBM and LAT), it is possible to study the γ-ray spectra of Fermi GRBs in unprecedented detail. In this work, we intend to conduct a comprehensive analysis of the most complete Fermi samples at present to systematically search for millilensing events. We will firstly apply the autocorrelation method to analyze light curves of each source in the sample, and make the preliminary screening of candidates in combination with hardness test (see Section 2). Then we will analyze the time-resolved spectrum for each candidate for further justification (see Section 3). Conclusion and discussion are presented at the end (see Section 4).

Data
In this work, we use the data observed by the Fermi GBM, which consists of 14 detector modules: 12 Sodium Iodide (NaI) detectors, covering the energies 8 keV -1 MeV, and two Bismuth Germanate (BGO) detectors, covering 200 keV to 40 MeV (Meegan et al. 2009). The data was downloaded from Fermi Science Support Center (FSSC) 's FTP site 3 , which contains 3000 GRBs up to April 2021. For each burst, we use the Time-Tagged Events (TTE) data for both spectral and temporal analysis, which records each photon's arrival time with 2 µs temporal resolution, as well as information regarding in which of the 128 energy channels the photon registered. The prebinned, eight energy channel (CTIME) data was used to carry out the background subtraction. In the temporal analysis (including the auto-correlation test and hardness test), we only use the data from the one NaI with the highest signal to noise ratio. In the spectral analysis, data from BGO detectors are involved.

Auto-correlation test
For each Fermi GRB, we produce its light curve from the TTE data with 0.01s resolution for short GRB and 0.1s resolution for long GRB. The background was determined by fitting first-order polynomials without the signal part. In order to search for two emission episodes with similar light-curve patterns, we apply the autocorrelation test to each GRB, with the following function: where σ is the standard deviation for the light curve f (t), ∆t is relative displacement for auto-correlation, and N is the bin number of the light curve.
For each burst, we can plot its a(∆t) − ∆t correlation curve. In principle, if there is no similar emission episodes in the light curve, the a(∆t) − ∆t curve should show a smoothly downward trend with the highest peak at ∆t = 0. Otherwise, if similar emission episodes indeed exist, multiple peaks would be superposed on the background decay component. Note that the residual noise may also lead to small wiggles in the a(∆t) − ∆t curve for dim GRBs, here we regard it as an effective peak only when the peak height is 3-σ higher than the background. 170 candidates, whose a(∆t) − ∆t curve contains at least one effective peak (peak at ∆t = 0 doesn't count), are selected as the preliminary candidates. ∆t of the highest effective peak could be taken as the possible millilensing time delay, whose uncertainty could be estimated by adding Poisson noise to the original lightcurve (Ukwatta et al. 2010;Hakkila et al. 2018;Veres et al. 2021, for a similar approach).

Hardness test
Hardness test has been widely used to justify the lensing effect (Paynter et al. 2021;Kalantari et al. 2021;Wang et al. 2021;Veres et al. 2021), based on the hypothesis that the flux ratio between gravitationally lensed pulses should not depend on energy (Paczynski 1987). In this work, we adopt three energy channels to carry out the hardness test: Low energy channel (8-50 keV), Medium energy channel (50-110 keV) and High energy channel (110-323 keV).
For each preliminary candidate, we first plot their light curves in different energy channels. Based on the auto-correlation results, we can divide each light curve into two similar episodes (we first selected the interval containing the first episode by visually identifying contiguous temporal bins with significant signal, and then selected the second interval with the same length as the first one, but with a certain time delay). Then we define hardness ratios HR HM and HR ML for each episode as where N is the total photon counts of the episode, B is the corresponding background photon counts, and i, j = {L, M, H} indicate the channel index. The uncertainty of the HR ij could be estimated as where Poisson noise is considered. Here, we require the preliminary candidate to pass the Hardness test only when their HR HM and HR ML for different episodes being consistent with the mean value within 1-σ region.

FINAL CANDIDATE SELECTION
In order to make a better justification on the final candidates, for each candidate, we have analyzed and compared the time integral spectrum and time-resolved spectrum of the two emission episodes in detail. The spectral fitting is performed by using Markov Chain Monte Carlo (MCMC) method with an automatic code "McSpecfit" , and in the analysis of time-resolved spectrum, the time bin for each slices are manually selected to ensure sufficient SNR, so as to give a reliable fitting. For each fitting, we adopt four different spectral models, including signal power-law, cutoff power-law, Band function and Blackbody function. We use the Bayesian Information Criteria (BIC; Schwarz (1978)) to test the goodness of each model, which shows that in most slices of all candidates, Band function and cutoff power-law models are clearly better than signal power-law and Blackbody models, while cutoff powerlaw performs slightly better than Band function. Therefore, we decided to use the best fitting parameters of cutoff power-law model to compare the similarity of the two emission episodes. We discuss each candidate in detail below.  ; trigger number 081126899) with a duration of 54.145 ± 0.923s (T 90 , the time interval between 5% and 95% of the cumulative flux). The analysis results are shown in Figure 1. Auto-correlation analysis suggests the time delay between two emission episodes is 30.6 ± 0.3. We select data from the T 0 − 1.7s to T 0 + 8.3s as the first episode and shift it to the second episode with ∆t = 30.6 s. The hardness ratios of two episodes are consistent both for the first episode (HR HM = 0.98 ± 0.10, HR ML = 0.59 ± 0.06) and for the second episode (HR HM = 0.84 ± 0.11, HR ML = 0.57 ± 0.07). Here we performed a so-called "χ 2 test" for testing the light curve similarity of two episodes, which considers the binned lightcurves of the two pulses as representing two distributions and asks if they are consistent with coming from the same parent distribution Mukherjee & Nemiroff (2021). Following the same algorithm as Mukherjee & Nemiroff (2021), we calculate the minimum χ 2 value for GRB 081126A, which is 100.94 for 100 degrees of freedom with the 0.1s time bin, corresponding to a p-value of 0.46.

GRB 081126A
The time-integrated spectrum parameters of two episodes are E peak1 = 336 +49 −26 , α peak1 = −0.69 +0.06 −0.10 and E peak2 = 238 +28 −21 , α peak2 = −0.61 +0.12 −0.10 . As shown in Figure 1, we divide each episode into 6 slides to compare the time-resolved spectrum. In most slides (4/6), the α and E peak uncertainties of two episodes both overlap within 1-σ region. At the same time, the spectrum parameters of two episodes are consistent within 2-σ confidence level in every slide for time-resolved spectrum. Considering that the detection of GRB 081126A has a relatively high SNR, and based on its excellent performance in light curve and spectrum analysis, we rank it as the first class candidate for millilensing event. Moreover, based on the early detection on the afterglow of GRB 081126A in the UV filters, it is suggested that GRB 081126A should have a redshift of approximately 2.8 < z < 3.8 Bhat & van der Horst 2008), which further increases the possibility for GRB 081126A being a millilensing event.
For a lensed GRB, one can estimate the redshift mass of the lens source with where f is the amplification ratio between two images, z L is redshift of the lens source, G is the gravitational constant, and c is the speed of light (Mao 1992). Here f could be calculated with the total counts and background counts of two episodes as f = (N 1 − B 1 )/(N 2 − B 2 ). For GRB 081126A, we have f = 1.35 ± 0.02. In this case, the redshift mass of the lens source can be estimated as M L (1 + z L ) = 5.1 +0.3 −0.3 × 10 6 M . Auto-correlation analysis suggests the time delay between two emission episodes is 42.1 ± 0.2s. We select data from the T 0 to T 0 +17s as the first episode and shift it to the second episode with ∆t = 42.1s. The hardness ratios of two episodes are consistent both for the first episode (HR HM = 0.63 ± 0.04, HR ML = 0.44 ± 0.02) and for the second episode (HR HM = 0.61 ± 0.07, HR ML = 0.44 ± 0.04). The χ 2 test obtains χ 2 min = 196.7 for 230 degrees of freedom with the 0.1s time bin, which corresponds to a p-value of 0.95 4 .
The time-integrated spectrum parameters of two episodes are E peak1 = 161 +8 −7 , α peak1 = −1.11 +0.03 −0.04 and E peak2 = 158 +13 −12 , α peak2 = −1.03 +0.06 −0.06 . As shown 4 It is worth noticing that the result of "χ 2 test" is very sensitive to the time bin sizes of the light curve. For instance, Mukherjee & Nemiroff (2021) found a relatively low p-value when the time bin for GRB 090717A is adopted as 1.024 s. Comprehensive analysis, especially the time-resolved spectral analysis, is thus essential for justifying millilensing effects. in Figure 2, for time-integrated spectrum, the α and E peak of two episodes are remarkably consistent within 1-σ confidence level. We then divide each episode into 10 slides to compare the time-resolved spectrum, and the spectrum parameters of two episodes are consistent within 1-σ confidence level in every slide except the α of slide 5, 6 and E peak of slide 2, 8, which are consistent within 2-σ confidence level. Considering its temporal and spectral analysis results, we rank GRB 090717A as the first class candidate that may have experienced millilensing effect, with the amplification ratio as f = 1.69 ± 0.01, and the redshift mass of the lens source as M L (1 + z L ) = 4.03 +0.05 −0.05 × 10 6 M .

GRB 081122A
GRB 081122A was triggered on 2008-11-22 12:28:12.211 UT (T 0 ) by Fermi-GBM (McBreen (2008); trigger number 081122520) with a duration of 23.30 ± 2.11s. The analysis results are shown in Figure  3. Auto-correlation analysis suggests the time delay between two emission episodes is 13.7±0.3s. We select data from the T 0 − 0.4s to T 0 + 4.0s as the first episode and shift it to the second episode with ∆t = 13.7s. The hardness ratios of two episodes are consistent both for the first episode (HR HM = 0.77 ± 0.06, HR ML = 0.67 ± 0.05) and for the second episode (HR HM = 0.63 ± 0.11, HR ML = 0.56 ± 0.07). The χ 2 test obtains χ 2 min = 97.98 for 88 degrees of freedom with the 0.05s time bin, which corresponds to a p-value of 0.21. Here 0.05s time bin is adopted because the total duration of GRB 081122A is relatively short. The time-integrated spectrum parameters of two episodes are E peak1 = 220 +16 −11 , α peak1 = −0.61 +0.06 −0.05 and E peak2 = 189 +37 −20 , α peak2 = −0.80 +0.10 −0.12 . As shown in Figure 3, we divide each episode into 6 slides to compare the time-resolved spectrum. The spectrum parameters of two episodes are consistent within 1-σ confidence level for each slide. The overall spectral analysis results provide positive evidence to support GRB 081122A to be a millilensing candidate. However, we notice that there are two peaks in the first emission episode but the light curve of the second episode seems no such a feature which may be due to the low SNR. We thus rank it as the second class candidate that may have experienced millilensing effect. With the the amplification ratio f = 2.22 ± 0.07, the redshift mass of the lensing source can be estimated as M L (1 + z L ) = 8.6 +0.4 −0.4 × 10 5 M .

GRB 110517B
GRB 110517B was triggered on 2011-05-17 13:44:47.600 UT (T 0 ) by Fermi-GBM (trigger number 110517573) with a duration of 23.04 ± 0.36s. The analysis results are shown in Figure 4. Auto-correlation analysis suggests the time delay between two emission episodes is 17.2 ± 0.1s. We select data from the T 0 − 4s to T 0 + 9s as the first episode and shift it to the second episode with ∆t = 17.2s. The hardness ratios of two episodes are consistent both for the first episode (HR HM = 0.51 ± 0.06, HR ML = 0.52 ± 0.05) and for the second episode (HR HM = 0.54 ± 0.07, HR ML = 0.48 ± 0.04). The χ 2 test obtains χ 2 min = 185.4 for 130 degrees of freedom with the 0.1s time bin, which corresponds to a p-value of 0.001.
The time-integrated spectrum parameters of two episodes are E peak1 = 112 +7 −5 , α peak1 = −0.45 +0.09 −0.09 and E peak2 = 103 +6 −5 , α peak2 = −0.43 +0.10 −0.10 . As shown in Figure 4, for time-integrated spectrum, the spectrum parameters of two episodes are consistent within 1-σ confidence level. After we divide each episode into 11 slides to compare the time-resolved spectrum, the α of two episodes are consistent within 1-σ region in every slide, but the E peak are only consistent within 1-σ region in 3 slides (27.3%), and the spectrum evolution of the E peak show different trend in two episodes , which cuts down the possibility for millilensing event. We notice that there are multiple peaks in both emission episodes, where the width distribution of these peaks is similar, but the peak amplitude distribution exists some differences in detail (that is why we get a small p-value in the χ 2 test). Considering its temporal and spectral analysis results, we rank GRB 110517B as the second class candidate for millilensing event. With the amplification ratio f = 1.04 ± 0.02, the redshift mass of the lensing source can be estimated as M L (1 + z L ) = 2.2 +1.1 −1.1 × 10 7 M .   It has long been proposed that GRBs have potential to be gravitationally lensed into multiple images, due to their high redshift nature. When the lens source has mass 10 4 ∼ 10 7 M (such as intermediate-mass black hole, compact dark matter, star cluster or population III star), the delay time between different images would be in order of 10 −1 ∼ 10 2 s, which is comparable with the typical duration of GRBs. In this case, signals from two images would be collected within one single trigger, manifesting as one GRB that contains two (or even more) emission episodes with similar properties in both temporal and spectral domain. Such kind of events are usually called millilensing events.
In previous works, a few GRBs have been proposed as the millilensing candidates, i.e., BATSE GRB 950830, Fermi GRB 200716C 5 and Fermi GRB 210812A. In this work, we conduct a systemic search for millilensing events from the Fermi GRB sample (3000 GRBs up to April 2021). We select preliminary candidates by performing auto-correlation test and hardness test, and then we analyze and compare the time integrated spectrum and time-resolved spectrum of the similar emission episodes. Eventually we find 4 interesting candidates which may have experienced millilensing effect, e.g. GRB 081126A, GRB 090717A, GRB 081122A and GRB 110517B. We rank GRB 081126A and GRB 090717A as the first class candidate based on their excellent performance both in temporal and spectrum analysis. We rank GRB 081122A and GRB 110517B as the second class candidates (suspected candidates), mainly because their two emission episodes show clear deviations in part of the time-resolved spectrum or in the time-integrated spectrum. The ratio of the number of lensed candidates to the total number of GRB in our sample is 2/3000 ∼ 4/3000. The redshift mass of the lensing source for our selected candidates are within the range of 8.6 × 10 5 M < M L (1 + z L ) < 2.2 × 10 7 M , which are very likely star clusters or super-massive black holes.
For a given GRB, the lensing optical depth can be written as (Zhou et al. 2021(Zhou et al. , 2022 where monochromatic mass distribution for lens source is assumed, f L represents the fraction of the mass of lens with mass M L to the total matter mass (dark matter mass + baryon mass), z S is the redshift of the GRB, z L is the redshift of the lens, n L is the comoving number density of the lens, H(z L ) is the Hubble parameter at z L , H 0 is the Hubble constant, and Ω m is the present density parameter of matter. Here we define the the maximum value of normalized impact parameter as y max (R f,max ) = R 1/4 f,max − R −1/4 f,max , by requiring that the flux ratio of two lensed images is smaller than a critical value R f,max (in this work we take R f,max = 5). If we accumulate a considerable number of GRBs with redshifts satisfying N (z S ) distribution, the integrated optical depth of all these GRBs would bē τ (f L ) = dz S τ (f L , z S )N (z S ).
Consequently, the expected lensing fraction for GRBs would be In this work, our results suggest that F = 2/3000 ∼ 4/3000 (assuming M L ∼ 10 6 M ), inferring f L = (0.5 ∼ 1) × 10 −2 . It is worth noting that although we have made a detailed analysis of each source in the sample, it can not be ruled out that there are still some possible candidates who have not been screened out. The main reason is that when the time delay is shorter than the intrinsic time scale of the burst, the signals of different images will overlap. This situation is not easy to be distinguished through auto-correlation analysis, especially when the signal itself has a complex structure. To be conservative, here we suggest to use 0.5×10 −2 as a lower bound estimate for the fraction of the mass of lens objects with mass M L ∼ 10 6 M to the total matter mass, inferring that the density parameter of lens objects with mass M L ∼ 10 6 M is larger than 1.5 × 10 −3 .
Up to now, all candidates (including our findings and previous work findings) are proposed solely based on the gamma-ray data analysis. However, the physical origin of gamma-ray burst (GRB) prompt emission is still poorly understood (Zhang 2018). It cannot be excluded that a GRB might have two emission episodes with inherently similar temporal and spectral characteristics (Lan et al. 2018). Multi-band observations would be essential to finally determine whether a GRB has really experienced the gravitational lensing effect (Chen et al. 2021). With the successful operation of many sky survey projects in multiple bands, such as all-sky gammaray monitors (e.g. Gravitational wave high-energy Electromagnetic Counterpart All-sky Monitor; ), sky survey detectors in the X-ray band (e.g. Einstein Probe; (Yuan et al. 2018)), and wide field of view monitoring system in optical band (e.g. Ground-based Wide Angle Camera system; (Wei et al. 2016)), more lensed GRBs are expected to be detected and accurately certified in the future.