Exploring Supermassive Compact Dark Matter with the Millilensing Effect of Gamma-Ray Bursts

In this work, we use the data released by the Fermi Gamma Ray Burst Monitor (GBM), which are downloaded from the Fermi Science Support Center’s (FSSC) FTP site ¹¹1https://heasarc.gsfc.nasa.gov/FTP/fermi/, and processed with the GBM data tools v1.1.1. 3000 GRBs were detected by Fermi/GBM up to August 1st, 2022. We use the time-tagged event (TTE) data for our following analysis. The TTE data file records each photon’s arrival time with 2 $\mu$s temporal resolution and which of the 128 energy channels the photon registered. In these currently availbale 3000 GRBs, a great deal of effort has been taken to search lensed GRB candidates. Here we collected their results as follow: Yang et al. [23] and Wang et al. [24] independently performed Bayesian analysis of the prompt emission light curves and energy spectra, concluded that GRB 200716C ($M_{\rm L,z}\approx 0.43\times 10^{6}\leavevmode\nobreak\ M_{\odot}$) is in favor of the millilensing scenario with two similar pulses. Meanwhile, Veres et al. [25] also carried out exhaustive temporal and spectral analysis for the sample and claimed that GRB 210812A (we use redshifted mass $M_{\rm L,z}\approx 1.13\times 10^{6}\leavevmode\nobreak\ M_{\odot}$ from N2 pulse model of GRB210812A, because the evidence using the N2 pulse model is more compelling than in the case of the N1 shape) shows strong evidence in favor of the millilensing effects. Later on, Lin et al. [26] obtained that four candidates pass the hardness test and showed similarities in both temporal and spectral domain, i.e. GRB 081126A ($M_{\rm L,z}\approx 5.1\times 10^{6}\leavevmode\nobreak\ M_{\odot}$), GRB 090717A ($M_{\rm L,z}\approx 4.02\times 10^{6}\leavevmode\nobreak\ M_{\odot}$), GRB 081122A ($M_{\rm L,z}\approx 0.86\times 10^{6}\leavevmode\nobreak\ M_{\odot}$), and GRB 110517B ($M_{\rm L,z}\approx 22\times 10^{6}\leavevmode\nobreak\ M_{\odot}$). GRB 081126A and GRB 090717A are ranked as the first-class candidates based on their excellent performance in both temporal and spectrum analysis [26]; 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 [26]. In our following analysis, we firstly assume that all these 6 candidates are millilensed GRB events for exploring properties of compact dark matter. However, it was argued that GRB 090717A, GRB 200716C, GRB 081122A, GRB 081126A, and GRB 110517A can not pass all millilensing tests (light curve similarity test and hardness similarity test) [28, 29]. Therefore, we also consider the case that only GRB 210812A in the 6 millilensing candidates is a real lensing system.

For the purpose of this work, we need to collect the distribution of redshift, signal-to-noise ratio (SNR), and the width of each main pulse in the whole GRB sample. In addition, we also need to collect the time delay and the lens mass for all lensing candidates. So firstly we conduct an initial screening of all the events to select GRBs with redshift (usually the redshift is obtained from the spectral lines of afterglow observations) and 171 GRBs are eventually selected. For these events, the time range of the main pulse is selected based on following criteria: 1)SNR of the pulse within the time interval should be greater than 30 and it is calculated as follow: ${\rm SNR}=(C_{\rm all}-C_{\rm bak})/C_{\rm bak}^{1/2}$, in which $C_{\rm all}$ is the photon counts in the time range, $C_{\rm bak}$ is the background photon counts in the time interval; 2)the gap of the pulse should be lower than its 10$\%$ height between other pulses. There are three kinds of situation for the selection. For the first situation, if GRB cannot satisfy the SNR or height criteria, we would choose the whole burst as the main pulse. For the second situation, if there is only one pulse of the GRB satisfying the SNR and height criteria, this pulse would be selected as the main pulse. For the third situation, if there are multiple pulses satisfying the SNR and height criteria, the earliest pulse would be selected as the main one. The width of the main pulse is defined similarly to $T_{90}$, i.e. time interval between $5\%$ and $95\%$ of the cumulative flux. As shown in Figure 1, we assume that the distribution of width and SNR of the main pulse for whole GRB sample is the same as the 171 GRB sample.

If we have the redshift information of these 6 millilensing candidates, we can write the redshift distribution of these lens in the form of the optical depth for each candidate as

$$P_{i}(z_{\rm l})=\frac{1}{\tau(z_{{\rm s},i})}\frac{d\tau(z_{{\rm s},i})}{dz_{\rm l}},$$

(1)

where $\tau(z_{{\rm s},i})$ is optical depth for each event and is expressed as

$$\begin{split}\tau(z_{{\rm s},i})=\int dm\int_{0}^{z_{{\rm s},i}}dz_{\rm l}\frac{dn(m,\Phi)}{dm}\times\\ \frac{d\chi(z_{\rm l})}{dz_{\rm l}}(1+z_{\rm l})^{2}\sigma(m,z_{\rm l},z_{{\rm s},i}),\end{split}$$

(2)

where $\chi(z_{\rm l})$, ${dn(m,\Phi)}/{dm}$, and $\sigma(m,z_{\rm l},z_{{\rm s},i})$ are the same as later Eq. (5) respectively. Unfortunately, we do not have redshift information for these millilensing candidates. As shown in the left panel of Fig. 2, we can only infer the redshift distribution of lens for 171 GRB sources by combining the redshift distribution of these GRB sources with Eq. (1). Then we combine this inferred redshift distribution of lens with the posterior distribution $p(M_{z,\rm L}|d_{i})$ for each millilensing candidate to infer $p(m|d_{i})$ as shown in the right of panel of Fig. 2. We find that the lens masses of 6 millilensing candidates lie in the range of $10^{5}-10^{7}\leavevmode\nobreak\ M_{\odot}$ ²²2We have tested that the redshift distribution of lenses has little influence on the estimation of magnitude of population information. Therefore, our methods and results should be reliable.. It is worth noting that here we use the skewnorm distribution to fit the posterior distribution $p(M_{z,\rm L}|d_{i})$ obtained from [26, 25, 24].

For population hyperparameters of compact dark matter $\Phi=[\bm{p}_{\rm mf},f_{\rm CO}]$ and $N_{\rm obs}$ detections of millilensed GRB events $d=[d_{1},...d_{N_{\rm obs}}]$, the likelihood follows a Poisson distribution without considering measurement uncertainty and selection effect

$$p(d|\Phi)\propto N(\Phi)^{N_{\rm obs}}e^{-N(\Phi)}.$$

(3)

However, with measurement uncertainty and selection effect taken into account, the likelihood for $N_{\rm obs}$ milllensed GRB observations can be characterized by the inhomogeneous Poisson process as [31, 30]

$$\begin{split}&p(d|\Phi)\propto N(\Phi)^{N_{\rm obs}}e^{-N_{\rm det}(\Phi)}\times\\ &\prod_{i}^{N_{\rm obs}}\int d\lambda L(d_{i}|m)p_{\rm pop}(m|\Phi),\end{split}$$

(4)

where the likelihood of one lensing event $L(d_{i}|m)$ is proportional to the posterior $p(m|d_{i})$. $N(\Phi)$ is the total number of millilensing events in the model characterized by the set of population parameters $\Phi$ as

$$\begin{split}N(\Phi)=\int dm\int dz_{\rm s}\int_{0}^{z_{\rm s}}dz_{\rm l}\frac{dn(m,\Phi)}{dm}\times\\ \frac{d\chi(z_{\rm l})}{dz_{\rm l}}(1+z_{\rm l})^{2}\sigma(m,z_{\rm l},z_{\rm s})N_{\rm s}P_{\rm s}(z_{\rm s}),\end{split}$$

(5)

where $\chi(z_{\rm l})$ is the comoving distance, ${dn(m,\Phi)}/{dm}$ is the comoving number density of the compact dark matter in a certain extended mass distribution $\psi(m,\bm{p}_{\rm mf})$

$$\frac{dn(m,\Phi)}{dm}=\frac{f_{\rm CO}\Omega_{\rm DM}\rho_{\rm c}}{m}\psi(m,\bm{p}_{\rm mf}),$$

(6)

where $\rho_{\rm c}$ is the critical density of the universe. In Eq. (5), $\sigma(m,z_{\rm l},z_{\rm s})$ is the lensing cross section

$$\sigma(m,z_{\rm l},z_{\rm s})=\frac{4\pi mD_{\rm l}D_{\rm ls}}{D_{\rm s}}y_{0}^{2}.$$

(7)

It should be emphasized that we take the impact parameter as $y_{0}=5$ because $y_{0}>5$ are difficult to be identified as lensing signals³³3If the SNR of the secondary peak is greater than 8, the SNR of the main peak is at least greater than 5000. Such a strong GRB signal is almost impossible.. In addition, $p_{\rm pop}(m|\Phi)$ is the normalized distribution of lens masses and written as

$$p_{\rm pop}(m|\Phi)=\frac{1}{N(\Phi)}\frac{dN(m,\Phi)}{dm}=\psi(m,\bm{p}_{\rm mf}).$$

(8)

Meanwhile, $N_{\rm det}(\Phi)$ is the number of detectable millilensing events and can be defined as

$$\begin{split}N_{\rm det}(\Phi)=\int dm\int dz_{\rm s}\int_{0}^{z_{\rm s}}dz_{\rm l}\frac{dn(m,\Phi)}{dm}\times\\ \frac{d\chi(z_{\rm l})}{dz_{\rm l}}(1+z_{\rm l})^{2}\sigma_{\rm det}(\lambda,m,z_{\rm l},z_{\rm s})N_{\rm s}P_{\rm s}(z_{\rm s}),\end{split}$$

(9)

where $\sigma_{\rm det}(\lambda,m,z_{\rm l},z_{\rm s})$ is the cross section that lensing signal can be detected and depends on the source parameters $\lambda=[{\rm SNR},w]$ via

$$\begin{split}\sigma_{\rm det}(\lambda,m,z_{\rm l},z_{\rm s})=\frac{4\pi mD_{\rm l}D_{\rm ls}}{D_{\rm s}}\times\\ [y_{\max}^{2}({\rm SNR})-y_{\min}^{2}(w)].\end{split}$$

(10)

The maximum value of the normalized impact parameter $y_{\max}$ can be determined by requiring that the two lensed images are greater than some reference value of flux ratio $R_{\rm f,max}$

$$y_{\rm max}(R_{\rm f,max})=R_{\rm f,max}^{1/4}-R_{\rm f,max}^{-1/4}.$$

(11)

We take the distribution of SNR in Figure 1 and reference value as function of SNR ($R_{\rm f,max}={\rm SNR}/8$) to obtain $y_{\max}$. In addition, the minimum value of the impact parameter $y_{\min}$ can be determined by the time delay of lensed signals $\Delta t$ and pulse widths $w$

$$\begin{split}\Delta t(m,z_{\rm l},y)=4m\big{(}1+z_{\rm l}\big{)}\times\\ \bigg{[}\frac{y}{2}\sqrt{y^{2}+4}+\ln\bigg{(}\frac{\sqrt{y^{2}+4}+y}{\sqrt{y^{2}+4}-y}\bigg{)}\bigg{]}\geq w.\end{split}$$

(12)

Since the minimum value of impact parameter $y_{\min}$ is not sensitive to the pulse width $w$ as shown in Eq. (12), we choose the average $\bar{T}_{90}\approx 27.6\leavevmode\nobreak\ \rm s$ of 171 GRBs to obtain $y_{\min}$ for our following analysis. Then the posterior distribution $p(\Phi|d)$ can be calculated from

$$p(\Phi|d)=\frac{p(d|\Phi)p(\Phi)}{Z_{\mathcal{M}}},$$

(13)

where $p(\Phi)$ is the prior distribution for population hyperparameters $\Phi$ and we assume a log-normal mass function for lenses

$$\begin{split}\psi(m,\bm{p}_{\rm mf}=[\sigma_{\rm c},m_{\rm c}])=\frac{1}{\sqrt{2\pi}\sigma_{\rm c}m}\times\\ \exp\bigg{(}-\frac{\ln^{2}(m/m_{\rm c})}{2\sigma_{\rm c}^{2}}\bigg{)}.\end{split}$$

(14)

Here, we set prior distributions for all the population hyperparameters $\Phi=[\sigma_{\rm c},m_{\rm c},f_{\rm CO}]$ as shown in Tab. 1. In addition, $Z_{\mathcal{M}}$ is both the normalized factor and the Bayesian evidence for the population model $\mathcal{M}$. This normalized factor can be calculated from the integral of the numerator of Eq. (13) over $\Phi$, i.e.

$$Z_{\mathcal{M}}=\int d\Phi p(d|\Phi)p(\Phi).$$

(15)