Fermi-GBM Observations of GRB 210812A: Signatures of a Million Solar Mass Gravitational Lens

Observing gravitationally lensed objects in the time domain is difficult, and well-observed time-varying sources are rare. Lensed gamma-ray bursts (GRBs) offer improved timing precision to this class of objects complementing observations of quasars and supernovae. The rate of lensed GRBs is highly uncertain, approximately 1 in 1000. The Gamma-ray Burst Monitor (GBM) onboard the Fermi Gamma-ray Space Telescope has observed more than 3000 GRBs making it an ideal instrument to uncover lensed bursts. Here we present observations of GRB 210812A showing two emission episodes, separated by 33.3 s, and with flux ratio of about 4.5. An exhaustive temporal and spectral analysis shows that the two emission episodes have the same pulse and spectral shape, which poses challenges to GRB models. We report multiple lines of evidence for a gravitational lens origin. In particular, modeling the lightcurve using nested sampling we uncover strong evidence in favor of the lensing scenario. Assuming a point mass lens, the mass of the lensing object is about 1 million solar masses. High-resolution radio imaging is needed for future lens candidates to derive tighter constraints.


INTRODUCTION
Strong gravitational lensing is a tool that serendipitously enhances our observing capabilities and offers new opportunities to study the Universe (see e.g. Oguri 2019). Gamma-ray bursts (GRBs) are energetic transient sources at cosmological distances, involving relativistic jets from stellar-mass black hole (BH) central engines. GRBs last from a fraction of a second to about 1000 s and typically show non-thermal spectra (see e.g. Kumar & Zhang 2015;Beloborodov & Mészáros 2017, for reviews). Given that the distance scale of GRBs spans a wide range (up to redshift z 9, Cucchiara et al. (2011)), a fraction of GRBs will show the imprints of strong gravitational lensing.
Strong gravitational lensing produces multiple images of the same source. The images differ in their intensity, but importantly their spectral shapes will be the same. Similarly, in time-varying sources, the temporal profile will be the same but shifted in time for different images (Schneider et al. 1992). The temporal and spectral invariance is the main defining feature of gravitational lensing. Lensed images are separated by up to arcseconds which are clearly below the resolution of current gamma-ray detectors. Gamma-ray instruments, however, have excellent time resolution, and temporal structures can be recorded with unparalleled accuracy (Meegan et al. 2009).
In one incarnation of the lensing scenario (Refsdal 1964;Rodney et al. 2021) applied to GRBs, also called macrolensing, a single event triggers the same instrument twice, and can be separated anywhere from a few hours to decades. The two triggers will have similar lightcurve shapes and spectra. The delay between the events scales linearly with the lens mass. Lens candidates for this scenario have masses in the 10 8 − 10 12 M range. Objects in this mass range include supermassive black holes (up to few ×10 10 M ), galaxies and clusters of galaxies. In practice, however, all traditional strong lensing time delay measurements are from galaxies or clusters of galaxies. GRB lensing rates are highly uncertain but somewhere on the order of 1 in 1000 GRBs should be affected (Mao 1992). After roughly 10,000 observed GRBs during three decades, no convincing macrolensing candidate GRB pair has been found (Nemiroff et al. 1994;Davidson et al. 2011;Veres et al. 2009; arXiv:2110.06065v1 [astro-ph.HE] 12 Oct 2021 Hurley et al. 2019;Ahlgren & Larsson 2020). The negative result is likely a combination of two effects: First, GRB detectors typically in low Earth orbit will miss a sizeable fraction of GRB lens echoes (see, however Hurley et al. 2019;Hui & MoonBEAM Team 2021, for existing and future all-sky instruments). Second, for weaker GRBs, with pulses close to the noise level, it is difficult to distinguish between the lensing scenario and just two unrelated but similar looking GRBs (Ahlgren & Larsson 2020).
In a different scenario, also called millilensing (Nemiroff et al. 2001, because the expected separation between the images is on the order of milli-arcseconds), the gravitational lens signature is imprinted upon the lightcurve of a single trigger. In this case, we have, e.g., two emission episodes with similar lightcurve patterns which can be separated by timescales spanning from a fraction of a second to a few minutes.
Recently, there has been an increase in claims of millilensing events. Paynter et al. (2021) presented convincing evidence for lensing in the short duration BATSE GRB 950830. Mukherjee & Nemiroff (2021a) raised some concerns based on the inconsistent flux ratio between the two pulses. Yang et al. (2021) and Wang et al. (2021) independently argued that the likely short GBM GRB 200716C shows milli-lensing signatures. Kalantari et al. (2021) reported on a different GBM lensing candidate, the long duration GRB 090717, selected based on the analysis of the auto-correlation function. The claim of Kalantari et al. (2021) was challenged by Mukherjee & Nemiroff (2021b) arguing that the two pulse shapes differ significantly. For the abovereported lensing candidates the flux of the first and the second emission episode is either at the same level or their ratio is 1.5.
We present observations of the long duration GRB 210812A and show that it is consistent with a gravitational lensing scenario. It is the first lensing claim with a flux ratio 3. We list multiple lines of evidence to support the lensing interpretation. We perform spectral and temporal analysis using Fermi-GBM data and complement it with additional data from INTEGRAL-SPI/ACS and Swift/BAT.
In Section 2, we present the observations, followed by tests for gravitational lensing in Section 3. We discuss the power of the tests in Section 5 and conclude in Section 6. Veres & Fermi-GBM Team 2021). Fermi-GBM consists of 12 NaI (referred to as n0, n1,..., n9, na and nb) and 2 BGO (referred to as b0 and b1) detectors covering the entire unocculted sky in the 8-1000 keV and ∼0.1-40 MeV energy range, respectively.
As reported by the automatic pipeline, the GRB location is RA= 39.7, dec= 69.7 degrees, with an error radius of 1.1 degrees (statistical only). At the trigger time, this position corresponds to a LAT boresight angle of 149 degrees. The GRB location was behind the spacecraft, meaning most GBM detector normals have a large angle to the source. Based on previous experience (e.g. Connaughton et al. 2016), this type of geometry results in a large number of detectors showing approximately equal count rates compared to ∼3 detectors with dominating signal for typical GRBs with small LAT boresight angle. Upon visual inspection of the NaI detectors' lightcurves in the 50-300 keV range, where GRBs are brightest, we find that all 12 NaI detectors detected the brighter first peak. Moreover, detectors n1 and n6 through nb also detected the fainter, second pulse. Both of the pulses were detected by the 2 BGO detectors.
Detectors n8, na, and nb showed the strongest signals. We used data from these detectors, along with b1, for spectral analysis. For the temporal analysis, we used detectors n1, n6 through nb, b0, and b1. For both spectral and temporal analysis, we used the 128 energy channel time-tagged event (TTE) data. We chose the pre-binned, 8 energy channel (ctime) data to carry out the flux-ratio test.
The targeted search (Blackburn et al. 2015;Goldstein et al. 2019) was designed to search for coherent sub-threshold signals (weak signals that did not trigger the instrument). As expected, we recovered both pulses of GRB 210812A with high significance. The search also provides a location for the burst (RA= 40.5 • , dec= 69.4 • ), consistent with the location of the automatic pipeline used for standard GRB analysis. We also find that the locations of the two pulses are consistent, meaning they do indeed belong to the same source.

Other observations
The location of GRB 210812A was outside the coded mask of Swift-BAT (Gehrels et al. 2004;Barthelmy et al. 2005) at the time of the trigger. The first peak of GRB 210812A is clearly present in the continuous 4channel data 1 , 2 , but the second peak is only discernible in the summed lightcurve.
INTEGRAL/SPI-ACS (Winkler et al. 2003;von Kienlin et al. 2003) detected GRB 210812A and clearly shows the two-peak structure. 3 We calculate the time delay between Fermi-GBM and ACS is due to the higher altitude of the INTEGRAL spacecraft to be dt ACS = −0.396 ms. We applied this correction to the ACS lightcurve in Figure 2.

Temporal properties
GRB 210812A consists of two pulses, separated by a quiescent period of about 30 s ( Figure 1). Using the GBM targeted search, we found no emission between the two pulses. For background estimates, we fit a third degree polynomial based on quiescent segments of the lightcurve before, after and in-between the two pulses. We also report no discernible emission around 33 s after the second pulse, indicating that this is not a periodic source. We further note that low level emission discernible in the summed lightcurve just before the start of the second pulse (T0+26.5 s to T0+29.6 s) is inconsistent with coming from the location of GRB 210812A, and thus it is unrelated. GRB 210812A has a duration of T 90 = 39.9 ± 3.6 s (10-1000 keV; T 90 marks the time interval between 5% and 95% of the cumulative flux) 4 . Analyzed separately, the two pulses have consistent duration: the duration of the first pulse is T 90,1 = 5.31 ± 0.68 s while the second  pulse is T 90,2 = 3.84 ± 1.64 s long. Taking the first pulse as a separate GRB, we classify it based on the GBM T 90 distribution (Bhat et al. 2016;von Kienlin et al. 2020), as a likely long GRB originating from the collapse of a massive star, with a probability of 87%. Conversely, the likelihood of GRB 210812A being a short GRB from a compact binary merger, based on the T 90,1 information, the observed distribution of Fermi-GBM GRBs and calculated consistently with Goldstein et al. (2017);Rouco Escorial et al. (2021) is 13 %.
Autocorrelation function -We measure the time delay between the two pulses using the autocorrelation function of the lightcurve between 10 and 1000 keV, with 64 ms resolution. To accurately determine the autocorrelation peak, we fit a 9th degree polynomial to the peak region. We estimate the uncertainty of the delay by adding Poisson noise to the original lightcurve (see e.g. Ukwatta et al. 2010;Hakkila et al. 2018, for a similar approach). After repeating the peak finding procedure, we take the uncertainty as the 1σ confidence region of the delays of the modified lightcurves and get: This is the first among the many time delay measurements between the pulses. We note here that this accuracy is typical of what one would expect for a lensed GRB separated by a longer timescale (macrolensing). The peak of the auto-correlation curve shows a 3.15σ excess over the smoothed curve calculated using the Savitzky-Golay filter. This excess satisfies the criteria of Paynter et al. (2021) for lensing, which classifies GRB 210812A as a lensing candidate for further scrutiny.
Spectral lag -The lag is a measure of the delay between high and low-energy photons (e.g. Norris et al. 2000). Typically, a GRB is harder initially, and the higher energy photons (100-300 keV) arrive earlier than the low energy (25-50 keV) photons. This relation is also used to estimate the redshift based on the empirical correlation between lag and luminosity (Norris et al. 2000). For the first pulse we find that the lag is τ lag,1 = 221.6 +123.4 −139.6 ms, while the second pulse has: τ lag,2 = −89.2 +255.9 −206.5 ms. The error on the second pulse is much larger, because of its weakness, especially in the 25-50 keV range.

Spectral Analysis
The spectrum of the first pulse is best fit by a power law with exponential cutoff, also known as the Comp-  tonized model (see Table 1 for the parameters and Figure 4 for the spectrum). The fluence in the 10-1000 keV range is The second pulse is weaker, and it is fit equally well by a simple power-law and the Comptonized model. We chose the Comptonized model, because it is more physical (the power law with photon index > −2 integrates to infinite energy). The difference in the goodness of fit measure (∆C-stat) when going from the simpler power law (2 parameters) to the Comptonized model (3 parameters) is ∼6. This is just below the ∆C-stat ≈ 8 used in e.g. Poolakkil et al. (2021) to select the more complex model. However, all the parameters of the Comptonized model are well constrained (Table 1), which justifies the use of this model. The fluence of pulse 2 is F 2 = (21.6 ± 2.6) × 10 −7 erg cm −2 .
We can also compare the time-resolved spectrum of the two pulses. Because of the relative weakness of the second pulse, we chose to fit the simple power-law model in bins of 0.256 s. We show the temporal evolution of the power-law indices for the two pulses and a linear fit to both as a function of time ( Figure 5). We shifted the second pulse by 33.3 s to highlight their similar behavior.

INDICATORS OF LENSING ORIGIN
Proving the gravitational lensing origin involves showing that the two pulses have identical pulse shapes, and their spectra are similar as well. It is conceivable that some physical mechanism produces two pulses with identical pulse shapes and spectra with no lensing involved. Without imaging the sources based solely on the gamma-ray observations, a lensing scenario is difficult to prove beyond the shadow of a doubt.  In this section, we will show, however, that the gamma-ray properties of GRB 210812A pass all the tests in the literature for a lensing origin. Furthermore, we apply the Bayesian evidence criterion that can select among models, and this method indicates strong evidence in favor of the lensing scenario.
The basic idea for establishing the statistical likelihood of lensing in the temporal domain is comparing two scenarios. First, we assume no lensing. We fit the pulse model and allow every parameter to vary freely. In the alternative scenario, where lensing is assumed, the second pulse is forced to have the same shape as the first one and differ only by a normalization factor and a time delay.
We use two types of pulses (Norris et al. 1996(Norris et al. , 2005. The first pulse with 5 parameters, referred to as "N1", is defined as: Here A is the amplitude at peak time, t max . σ r and σ d mark the rise and decay timescales of the pulse and ν is a shape parameter.
The second pulse with 4 parameters, "N2", has the following temporal dependence: where A is the amplitude, ξ is the asymmetry parameter, ∆ is the start time of the pulse, and τ is a duration parameter.

Indirect evidence
We explore a few properties of GRB 210812A that are not direct proofs of lensing, but they are necessary to any such claim.
Hard-to-soft evolution: Gamma-ray bursts typically become softer with time (e.g. Kaneko et al. 2006). This trend can be observed in the individual pulses of GRB 210812A as well (see Figure 5). For GRBs in general, the hard-to-soft evolution can be observed even across pulses. Specifically, for GRBs that show two pulses separated by a quiescent period, the second pulse is, in general, softer (e.g. Lan et al. 2018;Zhang et al. 2012). We find that the second pulse of GRB 210812A has the same spectral shape within errors or, subsequently, the same hardness, which is uncommon for typical GRBs.
Leading pulse is brighter: In case of simple lens models, the light ray traveling closer to the lens arrives later. It has a lower magnification than the first arriving light ray with the larger impact parameter (Krauss & Small 1991). We find that the first pulse is indeed visibly brighter and thus consistent with the expectation from a lensed source. We note that this criterion may not hold for complex lens models (Keeton & Moustakas 2009).

Spectrum of the pulses
For the spectral analysis we first selected the interval containing the second pulse by visually identifying contiguous temporal bins with significant signal. For the first pulse, we selected a source interval which is the same length as the second pulse (see Table 1 and Figure  2). The spectral fits of the two pulses yield consistent spectral shapes within errors: the peak of the energyper-decade or νF ν spectrum is E peak,1 = 324 ± 28 keV and E peak,1 = 283 ± 90 keV. To compare the two spectra, we plot the 1 and 2 σ confidence regions of the spectral shapes (Figure 4), accounting for the correlations between the parameters. We multiply the second pulse by the fiducial 3.5 number to show that the two spectral shapes are consistent.

Count ratio test
If the pulses are gravitationally lensed, the ratio between pulses should not depend on energy. Mukherjee & Nemiroff (2021a) investigated the count ratio of the two pulses as a function of energy for GRB 950830 and found a 2σ inconsistency. This test has the advantage that it is independent of the particular spectral shape.
For GBM we used the 8-channel ctime data to carry out this test on GRB 210812A. We considered all the detectors where the second pulse was visible in any channel (n1, n6 through nb, energy channels 1 through 6, and b0 and b1, energy channels 0 and 1), and data from ACS and BAT (25-350 keV). We find that the ratio of the pulses in all the channels and all three instruments is consistent with the mean value within 1.6 standard deviations (see Figure 5). We thus conclude that GRB 210812A passes the count ratio test for lensing.

χ 2 test
A simple and robust test that doesn't assume any pulse shape was introduced by Nemiroff et al. (2001). This test was recently applied to the claim of Kalantari et al. (2021) on GRB 090717 by Mukherjee & Nemiroff (2021b). Mukherjee & Nemiroff (2021b) conclude that based on the χ 2 test, the claim of gravitational lensing can be excluded at the 5σ level.
This test considers the binned lightcurves of the two pulses as representing two distributions and asks if they are consistent with coming from the same parentdistribution. After appropriately re-scaling the first pulse and taking into account the background, we perform a χ 2 -test for the hypothesis that the two lightcurves are drawn from the same distribution.
The test statistic is defined the following way: where the t is the time,r < 1 is the scaling factor, P i marks the background subtracted counts in the two pulses, and B i is the background counts (i = {1, 2}).
As an example (see Figure 3, top left), we consider data from T 0 − 2 s to T 0 + 7 s interval (first pulse) and compare it with the interval shifted by ∆t = 33.3 s (second pulse). We use 0.256 s resolution, summed over the detectors with good signal in the 45-300 keV range, and added the signal of the BGO detectors (120-400 keV).
We use the tte instead of ctime data, because the beginning of the first pulse has uneven temporal binning in the ctime data.
First, we find the minimum of the χ 2 expression in Equation 3 as a function ofr and we getr = 0.231 (see Figure 3, top-left). Next, we calculate the minimum χ 2 = 24.3 value for 34 degrees of freedom, which corresponds to a p-value of 0.89. Thus there is no statistically significant difference between the two distributions.
We compare the lightcurves of the two pulses for different temporal resolutions, energy ranges and instruments, in the panels of Figure 3. We consistently find there is no statistically significant difference between the two pulses. For example, Mukherjee & Nemiroff (2021c) performed a preliminary χ 2 analysis on GRB 210812A and concluded there is a ∼ 2.8σ discrepancy between the two pulses. Using the same energy range, and detector selection as Mukherjee & Nemiroff (2021c) (assuming in their notation detectors are numbered 1 through 12, and ctime energy channels 1 through 8) we performed the χ 2 test on 512 ms resolution lightcurves and in the energy range 22-800 keV (NaI detectors only, Figure 3, top-right). We findr = 0.231 minimizes χ 2 at a value of χ 2 = 11.5 for 16 degrees of freedom (p-value of 0.78), indicating the two lightcurves are consistent with being drawn from the same distribution. The BGO lightcurve (512 ms resolution) yields χ 2 = 14.027(dof = 16) and p-value of 0.597, and for ACS (0.4 s resolution) χ 2 = 24.662(dof = 22) and pvalue of 0.314 (bottom two panels of Figure 3).

Bayesian Model Comparison
Applying an idea from gravitational wave model selection Paynter et al. (2021) introduced the Bayesian evidence to compare the lensing scenario to the case where there is no lensing.
In the no-lens scenario we fit the lightcurve with two pulses, with all the parameters left to vary. We derive the Bayesian evidence Z NL by integrating the likelihood over the multi-dimensional parameter space. E.g. for the N1 pulse, this involves 10 parameters: I(t) = I N 1 (t|A 1 , σ r,1 , σ d,1 , t max,1 , ν 1 ) + I N 1 (t|A 2 , σ r,2 , σ d,2 , t max,2 , ν 2 ). In the lensing scenario the second pulse is constrained. It has the same shape parameters as the first pulse, only differing in the normalization (r) and the shift in the peak time (∆t), resulting in 5 + 2 parameters: I(t) = I N 1 (t|A, σ r , σ d , t max , ν)+I N 1 (t|A/r, σ r , σ d , t max +∆t, ν). Z L is the evidence in this case.
Formally, the Bayesian evidence (or simply evidence) has the following meaning (e.g. Speagle 2020): from Bayes' Rule, the probability of the model parameters Calculating the evidence Z is computationally intensive. It requires integrating the likelihood over a multidimensional parameter space. We use the bilby python package (Ashton et al. 2019) and dynesty nested sampler (Speagle 2020) to carry out the integration over the parameters to find Z. As in Paynter et al. (2021) the difference in ln Z is the natural logarithm of the Bayes Factor (ln BF = ln Z L − ln Z NL ) and it can be used to decide between the models. The Bayes Factor is additive, values from different, independent (e.g. different energy ranges) measurements can be added to perform model comparison.
We present the results of this analysis in Table 2. Conveniently, the nested sampling yields also the best fitting flux ratio and time delay. Data from every instrument and energy range where the second pulse was detectable provided positive Bayesian evidence in favor of the lensing scenario. The simplest mass model is the point mass lens, when the mass of the lens is concentrated in a projected region smaller than the Einstein radius of the source. In this scenario, we can derive the lens mass from the flux ratio and the time delay (e.g. Mao 1992): where z l is the lens redshift and M l is the lens mass, c is the speed of light, G is the gravitational constant.

MCMC lightcurve fitting
We fit the summed NaI and BGO lightcurves with the Markov-chain Monte Carlo method using the emcee (Foreman-Mackey et al. 2013) python package with both the N1 and N2 pulses (see Figure 1). This method cannot select between the lensing and no-lensing scenarios; however, it is fast and robust compared to the more computation-intensive nested sampling (see section 3.5). We can take the result of the lightcurve fits in the lensing scenario and derive the lens mass, M l (1 + z l ) in the point-mass approximation. The N1 pulse model leads to a flux ratio of r = 4.47 +1.06 −0.73 and delay time of ∆t l = 33.16 +0.19 −0.21 s. The corresponding point lens mass value is: (1 + z l )M l = 1.07 +0.16 −0.15 × 10 6 M .
For the N2 pulse the flux ratio is r = 4.19 +0.30 −0.25 and the time delay ∆t l = 33.11 ± 0.06 s. These lead to a point mass lens of mass: (1 + z l )M l = 1.13 ± 0.06 × 10 6 M .  Table 2. Bayes factor compilation for different instruments, energy ranges and pulse models using the bilby nested sampling method. The flux ratio and the time delay resulting from the fits are shown in the last two columns.

Singular Isothermal Sphere (SIS) lens model
This lens model is characterized by the line of sight velocity dispersion, σ v of its mass distribution. While in the point mass lens case, we could constrain the lens mass, here, it is only possible to restrict the velocity dispersion up to a distance scale D, where D = D OL D LS /D OS and D ij mark the angular diameter distance combinations between the observer (O), lens (L) and source (S), respectively. We assumed a GRB redshift z s = 1 and a lens redshift z l = 0.4.
A simple mass estimate based on the virial theorem yields a mass M ≈ 8 × 10 5 (σ v /15 km s −1 ) 2 (R/10 pc)M , where R = 10 pc is an assumed size considered typical for e.g. globular clusters for which the SIS model is a good approximation. The mass is broadly consistent with the point mass lens approximation.

DISCUSSION
In the previous section, we performed tests to confirm GRB 210812A is affected by strong gravitational lensing. Here we discuss the strengths of each test, analyze the unlensed GRB properties and present future detection prospects.

Spectrum
The most basic spectral test for a lensing scenario is the flux ratio test. Because gravitational lensing is achromatic, the flux ratio of the two pulses has to remain constant across energy ranges and different instruments. This is a robust measure because it doesn't depend on the assumed spectrum. The only caveat to consider if the GBM detectors' pointing has changed between the two pulses. In that case, the spectral responses change significantly between the pulses and the recorded counts cannot be compared across the emission episodes of GRB 210812A. The pointing of the detectors however has not changed by more than 5 degrees between the pulses, which means the detectors' response in the direction of GRB 210812A is essentially the same throughout the duration of the burst. The flux ratios across the energy ranges and instruments are clearly consistent with being equal (Figure 6). The weighted mean is 3.45, and we measure the most significant deviation for the ACS data point (green), which is 1.6 standard deviations away. Mukherjee & Nemiroff (2021c) showed a preliminary analysis of the hardness ratio (HR) for ctime channels 3 and 4 (4 and 5 in their notation) and claimed a 2.2 σ discrepancy between the HR of the two pulses. Formally, the count ratio (CR) test is equivalent to the HR test: the hardness ratio of Pulse 1 is HR(P1)=C(P1, Ch2)/C(P1, Ch1), where C(P1, Ch2) denotes the count rate in pulse 1 for channel 2 (higher energy channel for the hardness), and C(P1, Ch1), HR(P2) can be calculated analogously. The count ratio in the energy channel 1 is CR(Ch1) = C(P1,Ch1)/C(P2,Ch1). Thus from CR(Ch1)=CR(Ch2), it follows that HR(P1)=HR(P2). Because there are no strong outliers in the count ratio test, we expect the data to reflect this in the HR test. For the first pulse we find: HR= 1.405 ± 0.063 and for the second pulse: HR= 1.127 ± 0.155. We confirm the finding of Mukherjee & Nemiroff (2021c) that the second pulse indeed shows a lower HR. However, taking their difference and adding the errors in quadrature, we find that the discrepancy is only 1.66σ, which does not invalidate the lensing scenario.
Next, we carried out a spectral analysis of the two pulses using the GBM data. In the Swift-BAT data, the second pulse was only visible in the summed lightcurve, and INTEGRAL SPI-ACS had only data in 1 energy channel. Therefore we only used Fermi-GBM data for the detailed spectral analysis. The spectral parameters are consistent within errors (Table 1), and a plot of the spectral shapes (Figure 4) also shows that the two spectra overlap when considering the confidence regions. We thus conclude that the spectra of the pulses are consistent with the lensing interpretation. We note the advantage of the continuous 128 energy channel data of GBM over the 4 channel tte data available for BATSE, for which precise spectral fits were not feasible (Paynter et al. 2021).

Time history
Gamma-ray instruments have better temporal than spectral resolution. E.g., the number of spectral resolution elements in 10-1000 keV is 10, while the 5 s duration pulse with ∼0.1 s resolution has 50 temporal resolution elements, where 0.1 s is the approximate timescale of variations for GRB 210812A (Bhat et al. 2012). For this reason, the temporal study of the lensing scenario can provide more constraints.
The weaker second pulse is closer to the noise level of the detectors. Thus the shorter duration of the second pulse is in line with expectations for a weaker bursts with similar pulse shape. Nonetheless, the duration of the two pulses is still consistent within errors as expected from a lensing scenario.
Independent of any pulse models, we first performed a χ 2 test to compare the two lightcurves. The χ 2 test determines if the two lightcurves are consistent with being drawn from the same distribution. The χ 2 test showed that there is no significant difference between the lightcurves in different energy ranges and across instruments ( Figure 3). We note, however that weakness of the second pulse results in relatively large Poisson errors, and fine temporal structures in the second pulse, if there are any, are washed out. This somewhat reduces the power of this test, showing only that the general shapes of the pulses are consistent.
Next, we introduced pulse models from the literature and fit the lightcurve in different energy bands and different instruments. Using nested sampling, we evaluated the evidence in favor of the lensing scenario using the bilby code. Independent of the detector, energy range, or pulse model (Table 2), the evidence is consistently for the lensing scenario as opposed to the non-lensing scenario (BF>0 in all cases). We note that the logarithm of the Bayes Factor is not necessarily positive for the model with less parameters (the lensing model in our case). Indeed, e.g. Wang et al. (2021) (their Table 1) shows some cases with negative ln(BF).
The Bayes Factor differs depending on which pulse model we apply. We find that in energy ranges where the second pulse is relatively weak, the evidence for lensing is not as strong (but it still favors the lensing). The evidence using the N2 pulse model is more compelling than in the case of the N1 shape. This can be due to the larger number of parameters in the case of the N1 pulse shape.
We consider the NaI data alone, and the N2 pulse shape fits. The sum of ln (Bayes Factors) for the NaI detectors yields ln(BF) ≈ 10.7. Following Kass & Raftery (1995) and Thrane & Talbot (2019) we can assign colloquial meaning to this number. A difference of more than 8 is considered strong evidence in favor of the lensing model. We conclude that the Bayesian evidence thus supports the lensing interpretation. We note that the BGO and ACS lightcurves also provide additional evidence and, to a lesser extent, the Swift-BAT data as well.
The nested sampling provides both a selection criterion between the models through the Bayesian evidence and, at the same time, provides the parameters for time delay and flux ratio. We show the values in Table 2 and Figure 7. Different instruments, detectors, and energy ranges all yield a consistent solution, pointing to the lensing origin. Most importantly, we can derive the lens mass for both pulse shapes and arrive at a consistent picture indicating a ∼ 10 6 M lens (Figure 8). The MCMC method yields smaller errors than the nested samplig. This is due to the different fitting approaches of the two methods (see e.g. Speagle 2020, for details).

GRB properties
If GRB 210812A had not been lensed, its fluence would have been F 0 = F 1 1 − 1 r ≈ 6.2 × 10 −6 erg cm −2 , where F 1 is the fluence of the first pulse (see Section 3.2) and we took r = 4.5 for the numerical value. Using this flux value, we can get a broad range for the possible redshift of this GRB using empirical correlations between gamma-ray properties.
We scan the 0.1 to 5 redshift range and find that z 0.5 is consistent to within one sigma with the Amati relation (Amati et al. 2002) between the isotropic equivalent energy, E iso and the redshift-corrected peak energy, (1 + z)E peak . While this is a very crude estimate, it is in line with the average measured redshift of GRBs (Bagoly et al. 2006;Jakobsson et al. 2006) and it is consistent with our fiducial value of z s = 1.
The lag-luminosity relationship (Norris et al. 2000;Gehrels et al. 2006) provides another estimate of the redshift (L ∝ (τ lag /(1 + z)) −0.74 ). Taking the median value of the lag (221.6 ms) and scanning the 0.1 to 5 redshift range, we find the redshift of GRB 210812A within the range 0.9 < z < 1.5 is consistent at 1σ level with the lag-luminosity relation. While this is similarly an empirical relation, it further reinforces that a redshift of z s = 1 is a reasonable approximation.

Future events
Even though GRB 210812A had no multiwavelength follow-up, as more lensed GRB candidates are observed, we expect to have a well-localized counterpart eventually. Identifying the afterglow of a lensed GRB, showing two consistently fading images, would be the smoking gun evidence for the lensing scenario. The angular separation of the two lensed images on the sky are on the order of the Einstein radius, which for a 10 6 M point mass is θ E = ((4GM/c 2 )(D LS /D OL D OS )) 1/2 ≈ 3 mas(M/10 6 M ) 1/2 (D/0.6 Gpc) −1/2 (assuming z l = 0.4 and z s = 1, mas = "milli-arcsecond"). This falls just short of the resolution of 10m class optical telescopes (≈40 mas). The only conceivable way of resolving the two images is through very large baseline interferometry (VLBI) radio imaging (Casadio et al. 2021). The sensitivity of VLBI can be as low as 10 µJy (Venturi et al. 2020) and radio afterglows are detected in significant numbers at or above this flux (Chandra & Frail 2012). VLBI imaging of the two sources will provide additional information on the source and lens redshift and help constrain the lens model. Capturing a lensed GRB with a well-understood lens will allow fully exploit the accurate, sub-second time delay measurement achievable for GRBs and will allow for time-delay cosmography.

Lensing object
It is difficult to identify the type of lens that produced the two pulses of GRB 210812A. Possible objects include black holes or globular clusters. Populations of objects can be ruled out based on their number density and their contribution to the total lensing probability (Nemiroff 1989). The total number of millilensed GRBs can be estimated based on a search for lensing candidates in the entire Fermi-GBM GRB catalog. Our goal in this paper was to report only on GRB 210812A, and we leave population-level studies for a future work. Nonetheless, we can make a few general observations, based on e.g. the previously mentioned claims by Kalantari et al. (2021) and Wang et al. (2021); Yang et al. (2021). For 1-3 lens events and the total number of Fermi-GBM GRBs, N > 3100, the lensing rate is (3−9)×10 −4 . This is on par with the rate based on BATSE observations by Paynter et al. (2021).
A black hole mass of M ≈ 10 6 M lies at the lower end of the supermassive black hole population with measured masses (e.g. Woo et al. 2010) and at the upper end of the intermediate-mass black hole population (Greene et al. 2020). Without detailed counterpart observations it is unclear in which group, if any, the lens of GRB 210812A belongs to. Further lensed events however can provide essential constraints on the rates and origin of black holes in this mass range.
Globular clusters are similarly good lens candidates and their masses can indeed reach 10 6 M (Baumgardt & Hilker 2018). A SIS model is a good approximation to the velocity dispersion in a globular cluster. Paynter et al. (2021) found, however, that even globular clusters with one order of magnitude smaller mass, ∼ 10 5 M , do not exist in sufficient numbers to produce the 1/2700 rate of lensed GRBs for BATSE. In our case, we re-quire similar lensing probabilities but with ∼ 10 6 M globular cluster population. 10 6 M globular cluster lies above the approximately 2×10 5 M turnover mass in the mass function (Jordán et al. 2007). This means that the larger 10 6 M mass cannot compensate in total lensing probability the drop in number density. We thus conclude that the globular cluster lens can be tentatively ruled out, and a point mass lens e.g. a black hole is more likely. For more definitive statements on the nature of the lens precise localization and high resolution observations will be necessary.

CONCLUSION
In this paper, we presented multiple lines of evidence for GRB 210812A being gravitationally lensed. The two peaks in GRB 210812A have consistent spectrum, time profile, and spectral evolution. We determined the flux ratio and time delay with multiple methods and arrived at a consistent picture. The first pulse is approximately 4.5 times brighter and the delay between the pulses is 33.3 s. Assuming a point mass lens, this flux ratio and delay corresponds to a lens mass of (1+z l )M l = 10 6 M . There are only a few unchallenged claims in the literature for lensed GRB lightcurves. GRB 210812A presents the first strong evidence for lensing a long GRB with a flux ratio larger than 2. Future events will benefit from high resolution radio observations for definitive proof of lensing origin and detailed lens modeling.