Evidence for gravitational lensing of GRB 200716C

Observationally, there is a small fraction of Gamma-ray bursts (GRBs) with prompt emission observed by Fermi/GBM that are composed of two pulses. Occasionally, the distance to a GRB may be lensed when a high mass astrophysical object resides in the path between the GRB source and observer. In this paper, we describe GRB 200716C, which has a two-pulse emission and duration of a few seconds. We present a Bayesian analysis identifying gravitational lensing in both temporal and spectral properties, and calculate the time delay ($\Delta t\sim 1.92$ s) and magnification ($\gamma\sim 1.5$) between those two pulses based on the temporal fits. One can roughly estimate the lens mass to be about $2.4\times 10^{5}~M_{\odot}$ in the rest frame. We also calculate the false alarm probability for this detection to be about 0.07\% with trial factors, and a present-day number density of about $808 \rm~Mpc^{-3}$ with an energy density $\Omega\sim 1.4\times 10^{-3}$. If the first pulse of this GRB near the trigger time is indeed gravitationally echoed by a second pulse, GRB 200716C may be a short GRB candidate with extended emission.


Introduction
The theory of general relativity (GR) predicts that space is curved by compact objects, and the phenomenon arising from the deflection of electromagnetic radiation (light rays) toward the mass in a gravitational field is called gravitational lensing (Blandford & Narayan 1992). A point mass gravitational lens magnifies and makes two different images of the source when a massive object is located close to the line of sight between the observer and a source (see Treu 2010 for a review). The photons traveling a longer distance will arrive first, but those traversing a shorter path deeper into the gravitational potential of the lens will have a delayed arrival due to a larger time dilation. Thus, the gravitationally retarded image is dimmer than the first image (see Section 2 for details). The observational signature of such an effect is an initial pulse followed by a duplicate echoed pulse. The duration of the time delay between the two emissions depends on the mass of the gravitational lens and the magnification of the two images (Mao 1992;Paynter et al. 2021). The profile of the light curve of the two images should be similar even with their different intensities. However, the gravitational lensing process does not change the photon energies, such that all source images should have the same spectra (Paczynski 1987;Mao 1992).
Gamma-ray bursts (GRBs) are some of the most luminous and active high-energy transients that have been described since their discovery in 1963 (see Kumar & Zhang 2015 for a review), and their spectroscopically measured redshifts span a range from 0.0085 (Galama et al. 1998;Kulkarni et al. 1998) to 9.4 (Cucchiara et al. 2011) with more than 10 4 observed GRBs. The discovery of gravitationally lensed optical sources with redshifts ranging from 0.102 to 5.699, implies that GRBs may be gravitationally lensed occasionally (Paczynski 1986). If this is the case, GRBs play an important role in searching for evidence of gravitational lensing. Each image with a gravitationally induced time delay and different magnification can be detected through the observed burst light curve (Paczynski 1986;Blandford & Narayan 1992;Kalantari et al. 2021). Based on the time delay and the ratio of peak flux of the two images, one can roughly estimate the lens mass in the rest frame (Mao 1992;Paynter et al. 2021;Kalantari et al. 2021).
From an observational point of view, a small fraction of GRBs with prompt emission observed by the Fermi Gamma-ray Burst Monitor (GBM) are composed of two or more emission episodes with a quiescent time that may last up to ∼100 s in the rest frame (Koshut et al. 1995;Lazzati 2005;Burlon et al. 2008;Bernardini et al. 2013;Hu et al. 2014;Lan et al. 2018). More interestingly, Lan et al (2018) performed a systematic analysis of both the spectral and temporal properties of GRBs with prompt emission observed by Fermi/GBM showing two-episode emission components in the light curves with quiescent times of up to hundreds of seconds. Statistically speaking they found that the spectral properties of the two-episode emission components are not significantly different, but they did not analyze carefully the light curves of those two-episode components. Recently, Paynter et al. (2021) claimed that they have found a possible signature of a gravitational lens in the light curve of GRB 950830 with two-episode emission. This could mean that the two-episode emission signatures are gravitationally lensed images of the same single-episode source. However, they do not present more details of the spectral properties of the two-episode emission.
One question is whether we can search for robust signatures of gravitational lensing in GRBs that produce two images within the source-lens-observer geometry and manifest in both the light curves and spectra. By systematically searching for more than 3000 GRBs observed with both Fermi/GBM and the Swift Burst Alert Telescope (BAT), we found GRB 200716C with two-episode emission. Its temporal and spectral properties satisfy the requirements of the theoretical predictions of gravitational lensing. In this paper, we show the basic theory of gravitational lensing in §2. Then, we discuss the evidence for gravitational lensing of GRB 200716C based on the observational data. A comprehensive data reduction and analysis of GRB 200716C is presented in §3, and a lens mass estimate is shown in §4. Conclusions are drawn in §5 with some additional discussion.

Basic theory of gravitational lensing
Considering a light ray from a distant source approaching a point mass (M), the bend angle α in the geometric optics limit is given as where b is the impact parameter by denoting the distance of closest approach of the ray to the mass, and G and c are the gravitational constant and speed of light, respectively. Figure  1 is a cartoon picture of the point mass gravitational lens geometry. First, let us label the observer-source distance D os , the observer-lens distance D ol , and the lens-source distance D ls . By assuming the weak field and thin-lens approximation, one has α = 4GM c 2 b ≪ 1 (i.e., for a small angle) and b ≪ D ol which implies θ ≪ 1. Since β < θ, β is also small angle. Based on the small angle mathematical geometry of projecting on a vertical line, we can write Combining Eq. (1) and Eq.
(2), one can solve the quadratic equation for θ, and find two solutions, For small angles, one multiplies both sides of Eq. (3) by D ol to solve for b, Thus, there will always be two lensed images for a point mass lens (Blandford & Kochanek 1987).
In order to find out the relationship between time delay (∆t) and magnification (γ) from the unlensed to the lensed system, we define a critical radius (also called the Einstein radius), inside which significant magnification takes place because the lensing changes the cross section but not the surface brightness (Turner et al. 1984). By defining a dimensionless impact parameter f = λ/r cr , the Eq.(4) can become The magnification (or the ratio of fluxes of individual images) can be expressed as The time delay is contributed by two effects when the arrival of photons follows the two paths in Figure 1. One is geometric due to different path lengths and the other is that two rays experience different general relativistic time dilations when the two paths traverse different gravitational potentials (Weinberg 1972). Thus, the time delay can be given as where M z = M(1 + z) is the redshifted lens mass.
GRBs have a good temporal resolution in the γ-ray band, and the time delay and magnification between the two images can be observed by considering both the difference in geometric path and the relative difference in gravitational field strength. Thus it is easy to estimate the mass of the gravitational lens:

Data reduction and analysis
In order to test how many of the GRBs observed by Fermi/GBM are potentially gravitationally lensed, as of 2021 July, we downloaded the original GBM data (12 NaI and 2 BGO detectors) of 3035 GRBs from the public science support center at the official Fermi website 1 . We employ the Bayesian Block algorithm to identify the light curves, and extract the spectrum using our automatic code "McSpecfit". Please refer to our previous paper (Lan et al. 2018) for more details on data analysis with the Bayesian Block algorithm, and to Zhang et al. (2018) for details on the spectral fitting. There are two criteria adopted for our sample selection. First, the GRB prompt emission must have two-episode (or more) emission, and the signal-to-noise ratio (S/N) of the emission episodes should be greater than 3σ. Second, the spectra of the two-episode (or more) emission should be similar to each other. After searching 3035 GRBs, we find that only GRB 200716C satisfies our criteria.

The basic observations of GRB 200716C
GRB 200716C triggered Swift/BAT, Insight-HXMT, and Fermi/GBM. Due to the lack of public Insight-HXMT data, in this section we focus on introducing the prompt emission of GRB 200716C observed by Swift/BAT and Fermi/GBM, as well as the afterglow (both X-ray and optical) observations. GRB 200716C triggered the BAT at 22:57:41 UT on 16 July 2020 (Ukwatta et al. 2020). We downloaded the BAT data from the Swift website 2 , and use the standard HEASOFT tools (version 6.28) to process the BAT data. For more details of the analysis, please refer to Sakamoto et al. (2008); Zhang et al. (2009);and Lü et al. (2020). The light curves in different energy bands are extracted with the time-bin size 8 ms. Then, we calculate the cumulative distribution of the source counts using the arrival time. The light curve shows two prominent peaks with a duration of about 5.3 s in 15-150 keV (see Figure 2), but weak activity is still visible until about 90 seconds.
At 22:57:41.18 UT on 2020 July 16, the GBM was triggered and located GRB 20716C (Veres et al. 2020). GBM has 12 sodium iodide (NaI) and two bismuth germanate (BGO) scintillation detectors covering the energy range from 8 keV to 40 MeV (Meegan et al. 2009). We downloaded the corresponding Time-Tagged-Event data from the public data site of Fermi/GBM 3 . For more details of the light-curve data reduction procedure refer to Zhang et al. (2016). The light curves of the n0 and b0 detectors with 8 ms and 64 ms time bins are shown in Figure 3, and consist of two pulses with a duration 3.3 s in 50-300 keV. There is no significant weak emission after the second pulse in the GBM temporal analysis.
The X-ray telescope (XRT) began observing the field at 22:59:04.2 UT, 82.9 s after the BAT trigger (Ukwatta et al. 2020). We made use of the public data from the Swift archive 4 (Evans et al. 2009). The X-ray light curve seems to be a power-law decay until ∼ 10 5 s with decay slope α 0 = 1.55 ± 0.02 (see Figure 4). Kann et al. (2020) observed the position of the afterglow with the 1.23 m Calar Alto telescope starting with the second Swift orbit and found that the decays follow a broken power law with decay slopes α 1 = 0.8 ± 0.04, α 2 = 5.5 ± 1.3, and break time t b = (3.8 ± 0.26) × 10 4 s (see Figure 4).

Light-curve fits of GRB 200716C
The light curve of GRB prompt emission with pulses is usually described with the fastrise exponential-decay (FRED) model (Norris et al. 1996). In order to test the consistency of structure for the two pulses, we also employ the FRED model to fit the pulses of GRB 200716C. By invoking the public code from Paynter et al. (2021), we used the same method from Paynter et al. (2021) to fit the light curve 5 . Here, we adopt two approaches to fit the data. Firstly, we used the same parameters (except the peak time and normalization) of one FRED model to fit the two pulses and obtain the values ln(Z L ), if we believe they are gravitationally lensed (called "FL"). Next, we used two FRED models to fit the two pulses with different parameters to get the values ln(Z NL ), if they are independent of one another (called "FF"). The light curve of a statistically significant gravitational lensing candidate GRB 200716C is shown in Figure 2. The reconstructed curves of the best model fit are plotted in black. We also present the difference between the true light curve and the posterior predictive curve in different energy channels. We find that the residuals are consistent with zero, which means the lens model we selected is a good one. On the other hand, in order to determine which model is preferred by the data we also calculate the Bayesian evidence for each model with the Bayes factor (lnBF ), which is defined as ln(BF ) = (ln(Z L ) − ln(Z NL )). A ln(BF ) that is larger than 8 is considered strong evidence for supporting one model over another (Thrane & Talbot 2019;Paynter et al. 2021).
We separate the Swift/BAT light curves into four available broadband energy channels, and independently calculate the value of ln(BF ) in those four channels (see Table 1). We find that the values of ln(BF ) are between −0.1 and 7.0 in each channel, and the total ln(BF ) value from each of the channels is about 15.24 in favour of the lensing hypothesis. This is strong statistical evidence supporting the lensing hypothesis.
Similar to the pulse fitting of Swift/BAT data, we also apply the FRED model to fit the the Fermi/GBM data. The reconstructed curves of the best model fits are plotted in black (see Figure 3). The residual in the different energy channels are also consistent with zero, indicating that the lens model is the best one. Here, we calculate the Bayes factor in four available energy channels (see Section 4) with 8 ms and 64 ms time bins, respectively. For the 8 ms time bin, the values of ln(BF ) are between 0.5 and 9.0 in each channel (see Table 1), and the total ln(BF ) value from each of the channels is about 19.94 in favor of the lensing hypothesis. But for the 64 ms time bin, the ln(BF ) is −0.5 during the first energy channel and in the other three channels it ranges from 4.0 to 9.0. The total ln(BF ) value from each of the channels is about 19.56, which is close to the value of ln(BF ) for the 8 ms time bin. This suggests that the total ln(BF ) value for each energy channel seems to be not dependent on the time resolution. At the least this is also strong statistical evidence supporting the lensing hypothesis.

Extracting and fitting the spectrum of GRB 200716C
We do not extract the spectrum of GRB 200716C observed by BAT due to its narrow energy band, but focus on the wide energy band in GBM. We extract the time-averaged spectrum of the first (time interval (−0.3 − 1.9) s) and second (time interval (1.9 − 4.1) s) pulses of GRB 200716C, respectively. The background spectra are extracted from the time intervals before and after those two pulses and modeled with an empirical function (Zhang et al. 2011). The spectral fitting is performed by using a Markov Chain Monte Carlo (MCMC) method with our automatic code "McSpecfit" in Zhang et al. (2018). We adopted several spectral models, which we usually select to test the spectral fitting of a burst, i.e., power law (PL), cutoff power law (CPL), Band function (Band), and Blackbody (BB), as well as combinations of any two models. Then, we compare the goodness of the fits of the two pulses, respectively (see Table 2). Invoking the Bayesian information criteria (BIC; Lü et al. 2017), we find that the CPL model is the best model that adequately describes the observed data. The CPL model fit is shown in Figure 5, as well as the parameter constraints of the fit. For the first pulse, it gives peak energy E p,1 = (524 ± 97) keV, and a lower energy spectral index of α 1 = 0.96 ± 0.05. For the second pulse, one has E p,2 = (566 ± 164) keV, and α 2 = 0.98 ± 0.08. The best-fit parameters of the CPL fits and other models are listed in Table 2.
Within the error range in the spectral data, the spectral properties of the two pulses in the CPL model are consistent with one another. This consistency is a the prediction of the lensing hypothesis. Based on the above analysis, both the light curve and spectral properties support that GRB 200716C is gravitationally lensed.

Estimating the lens mass of GRB 200716C
In order to determine whether the two pulses of GRB 200716C are a false alarm, based on the method of Paynter et al. (2021), we also calculate the false alarm probability 1 − p lens = 1 1+ln(BF)/N , where N = 3035 is total number of GRBs observed by Fermi/GBM. One has 1 − p lens = 7.3 × 10 −4 with the 8 ms time bin. In other words, the false alarm probability for this detection is about 0.07% with trial factors. Moreover, we also calculate that number density is about 808 Mpc −3 with an energy density Ω ∼ 1.4 × 10 −3 by assuming a redshift for GRB 200716C of z = 0.348 (D'Avanzo & CIBO Collaboration. 2020) and the average redshift of GRBs observed by Swift z ∼ 2.2 (Xiao & Schaefer 2011).
The gravitational lens will not change the photon energies when the photons travel close to compact objects, which means that all wavelengths of the light curve are equally affected by gravitational fields. In other words the time delay of different pulses is independent of the photon energy and it should be the same in different energy channels. Also, the gravitational magnification of each image is identical for every wavelength. In order to test this hypothesis with the observed data, we separate the Swift/BAT and Fermi/GBM light curves into four available broadband energy channels, respectively 6 .
Based on the light-curve fits for each energy channel and adopting a method similar to Paynter et al. (2021), one can easily to calculate the time delay and magnification. For Swift/BAT data, we roughly calculate ∆t ∼ 1.93 s and γ ∼ 1.54. For the 8 ms time bin of Fermi/GBM data, one has ∆t ∼ 1.92 s and γ ∼ 1.49. For the 64 ms time resolution, one has ∆t ∼ 1.92 s and γ ∼ 1.52. This indicates that both time delay and magnification are also independent of the time resolution. Figure 6 shows the peak flux ratio as a function of energy channels for prompt emission observed by Swift/BAT and Fermi/GBM (8 ms and 64 ms time bin). The ratio seems to be consistent across the different energy channels and time bins. By invoking the Eq. (10), as well as adopting ∆t ∼ 1.92 s and γ ∼ 1.5, one can roughly estimate that the lens mass in the rest frame is about 2.4 × 10 5 M ⊙ . There are several astrophysical objects within this mass range, such as globular clusters, diffuse galaxies, dark matter, and black holes (Paynter et al. 2021).
If GRB 200716C was lensed by a globular cluster, the estimated cosmic energy density of globular clusters Ω gc ∼ 8 × 10 −6 should be consistent with that of predictions. However, it is inconsistent as we infer energy densities much larger than that of globular clusters (see Figure 7). If the astrophysical object is a diffuse galaxy then it should have strong γ-ray and radio emission, which is inconsistent with current observations (Mihos et al. 2005). The other possible astrophysical object is an intermediate-mass black hole (Paynter et al. 2021), but whether black holes in this mass range exist remains an open question. By comparison with the result of Paynter et al. (2021), we find that the inferred lens mass of GRB 200716C is about 4 times higher than that of GRB 950803, and the inferred energy density of GRB 200716C is also about 3 times larger than that of GRB 950803. This result is consistent with that of Paynter et al. (2021).

Conclusion and discussion
GRB 200716C was observed by Swift, Fermi, and Insight-HXMT to have a duration of few seconds. The prompt emission of this GRB consists of two pulses and weak emission (called "extended emission") lasting ∼ 90 s after the second pulse is visible in the Swift/BAT, but not visible in the Fermi/GBM temporal analysis. In this paper, we presented a comprehensive analysis of its temporal and spectral data, and tested whether the first pulse of GRB 200716C near the trigger time is indeed gravitationally echoed by a second pulse, indicating that both pulses are gravitationally lensed images of the same single source pulse.
Firstly, we separated the Swift/BAT and Fermi/GBM light curves into four available broadband energy channels, respectively. The FRED model is invoked to fit the profile of two pulses in each channel by adopting the public code from Paynter et al. (2021). In Figure 2, we plot the light curve of 200716C observed by Swift/BAT, as well as the FRED model fits in different energy channels. The model fit is used to subtract from the true observed light curve and obtain the residuals. We find that the residuals are consistent with zero, which means that the lens model we selected is favored. Then, we independently calculate the Bayesian evidence for each model with Bayes factor (ln(BF )). We find that the total ln(BF ) value from each of the channels is about 19 for BAT and GBM (even with different time resolution). This value is much larger than 8, and so the lensing hypothesis is favored. It is also independent of the time resolution of the prompt emission. Moreover, we also extract the spectral parameters by using the MCMC method with our automatic code "McSpecfit" in Zhang et al. (2018). Several spectral models (PL, CPL, Band, and BB), or even combinations of any two models, are selected to fit. We find that the CPL model is the best one that adequately describes the observed data by comparing the goodness of the fits of the two pulses, respectively. Both the E p and α values of those two pulses are consistent with one another within the error range. This consistency is a prediction of the lensing hypothesis and is strong statistical evidence to support for the lensing hypothesis of GRB 200716C.
One basic question is whether the lensing signal from GRB 200716C is a false alarm. In order to test this question, we calculate the false alarm probability for this detection, which is about 0.07% with trial factors based on the method of Paynter et al. (2021). By adopting the redshift of GRB 200716C to be z = 0.348 and the average redshift of GRBs observed by Swift to be z ∼ 2.2, we estimated the number density as 808 Mpc −3 with an energy density Ω ∼ 1.4 × 10 −3 . On the other hand, we adopted a method similar to Paynter et al. (2021) and after making light-curve fits for each energy channel, we calculated the time delay and magnification of the pulses to be ∆t ∼ 1.92 s and γ ∼ 1.5, respectively. We find that the time delay and magnification of the two pulses are independent of the time resolution of the light curve. The inferred lens mass is about 2.4 × 10 5 M ⊙ , which is a mass consistent with several astrophysical objects such as globular clusters, diffuse galaxies, dark matter, and black holes (Paynter et al. 2021). However, the globular clusters and diffuse galaxies seem unlikely to be the candidate astrophysical objects. The black hole is a potential candidate, but more observations are needed to confirm this in the future.
Upon finishing this paper, our attention was drawn to Wang et al. (2021), who performed an independent analysis on GRB 200716C to discuss the same points. We find that there are two points of difference between this paper and Wang et al. (2021). First, the spectral fitting results of the two pulses are different, which may be caused by the different time interval selected and different fitting methods for the two papers. Several spectral models (PL, CPL, Band, and BB), or even combinations of any two models, are selected as fitting functions in our paper by using the MCMC method in our automatic code "Mc-Specfit". Wang et al. (2021) used only the Band function and CPL models to do fits but did not invoke an MCMC method to do that. Second, the estimated lens mass is slightly different for the two papers, but within the same order. The reason for this may be the selection of different time delay and magnification values. We used the average time delay and magnification values of Fermi and Swift in different energy bands to roughly estimate the lens mass, but Wang et al. (2021) presented the time delay and magnification values in each energy band and then estimated the lens mass.
If the GRB 200716C is indeed gravitationally lensed, the total duration of the prompt emission of this GRB should be the duration of any one pulse. If this is the case, then GRB 200716C should be a typical short-duration GRB with extended emission. Wang et al. (2021) claim that the E p − E γ,iso of GRB 200716C is located in the population of typical short GRBs, even for individual pulses by assuming a possible 7 redshift z = 0.348 (D'Avanzo & CIBO Collaboration 2020). At least for this case, due to the lack of accurate information on emission or absorption lines in the spectrum, we only can only find some indirect evidence for the gravitational lensing of GRB 200716C. The "Smoking gun" of gravitational lensing of GRBs is not only the consistency of the temporal and spectral properties with predictions from gravitational lensing, but the consistency with some empirical relations, and indeed accurate information of its host galaxy with two images. With the improvement of detection technology, we encourage observers in the future to invoke large optical telescopes to followup, especially for the GRBs with two-pulse emission. Moreover the light-curve behaviors between the X-ray and optical are quite different, so it makes this an event of interest. We also need to carry out a follow-up in the future.
Since the lensing signal could be due to similar-looking pulses of the GRB, the lensing hypothesis is one possible explanation for the double-pulse structure of GRB 200716C. On the other hand, the double pulse associated with a GRB 200716C-like event or even repeating pulses could be an intrinsic feature of the GRB prompt emission (Veres et al. 2021). In this case, it would be impossible to confidently detect lensing by looking at the similarity of the pulses.
We acknowledge the use of the public data from the Swift data and Fermi data archive. This work is supported by the National Natural Science Foundation of China (grant Nos. 11922301, and 12133003), the Guangxi Science Foundation (grant Nos. 2017GXNSFFA198008, and AD17129006), the Program of Bagui Young Scholars Program (LHJ), and special funding for Guangxi distinguished professors (Bagui Yingcai and Bagui Xuezhe).    . The bottom four panels correspond to residuals that show the data after the template has been subtracted for different energy channels. The colored shaded regions are the 1σ standard statistical error. These panels seem to show that the lens model is a reasonable fit.