The Luminosity Distribution of Short Gamma-Ray Bursts under a Structured Jet Scenario

The joint detection of gravitational wave (GW) and electromagnetic radiation from the binary neutron star merger event GW170817 marks a breakthrough in the field of multi-messenger astronomy. The short gamma-ray burst (sGRB) GRB 170817A, associated with this binary neutron star merger event, has an isotropic-equivalent gamma-ray radiation luminosity of 1.6 × 1047 erg s−1, which is much lower than that of other sGRBs. The measurement of the superluminal movement of the radio afterglow emission confirms the presence of the relativistic jet, and the emission features can be well explained by the structured jet model. In this paper, we calculate the luminosity distribution of sGRBs and its evolution with redshift based on the structured (Gaussian) jet model, and find that the typical luminosity increase with redshift, for nearby sGRBs (such as for luminosity distance less than 200 Mpc) the typical gamma-ray luminosity is just around 1047–1048 erg s−1, which naturally explains the very low radiation luminosity of GRB 170817A. We derived the detection probability of sGRBs by Fermi-GBM and found that the expected detection rate of sGRBs is only about 1 yr−1 within the distance of several hundred Mpc. We explored the effect of the power-law index α of the merger time distribution on the observed characteristics and found that it had little effect on the observed luminosity and viewing-angle distributions. However, it is very interesting that, for different values of α, the distributions of the number of observed sGRBs are quite different, so it is possible to determine the value of α through observed distributions of the number of sGRBs. We used the Bayesian method to make a quantitative analysis and found that the value of α may be identified when the number of observed sGRBs with known redshifts is more than 200. Finally, we compare our results of gamma-ray luminosity distribution with sGRBs with known redshifts, and found that our results are consistent with the observation, which implies that our simulation results can reproduce the observed luminosity distribution well.


Introduction
Gamma-ray bursts are the most violent explosions in the universe. Based on their durations, gamma-ray bursts are usually classified into two types, long gamma-ray bursts (LGRBs) and short gamma-ray bursts (sGRBs), which have different origins (for reviews, see Mészáros 2006;Zhang 2007;Gehrels et al. 2009). LGRBs are usually thought to originate from the collapse of massive stars (Woosley 1993). For a long time, sGRBs were only associated by indirect evidence with the merging of double compact stars (Nakar 2007); the indirect evidence included the location of sGRBs in host galaxies, nondetection of associated supernova (SN), large galaxy offsets, weak spatial correlation between sGRBs and star formation regions within their host galaxies (Nakar 2007;Berger 2014, for reviews), and, in particular, the identification of the socalled Li-Paczynski macronovae/kilonova in GRB 130603B (Berger et al. 2013;Tanvir et al. 2013), GRB 060614 Yang et al. 2015), GRB 050709 (Jin et al. 2016) and recently in GRB 070809 (Jin et al. 2020). The joint detection of gravitational wave (GW) and gamma-ray emission from the binary neutron star merger event GW170817/GRB 170817A marks a breakthrough in the field of multi-messenger astronomy, and provides the most conclusive evidence that at least a fraction of sGRBs are indeed from the merger of double neutron stars (Abbott et al. 2017a(Abbott et al. , 2017bGoldstein et al. 2017;Savchenko et al. 2017).
However, for GRB 170817A, the first sGRB associated with a GW signal, its prompt gamma-ray and afterglow emission show some peculiar characteristics, which are different from other sGRBs, for example, the isotropic-equivalent prompt gamma-ray emission energy is just about 3×10 46 erg and its luminosity is just about 1.6×10 47 erg s −1 , much lower than that of typical sGRBs. In addition, the X-ray and radio afterglows first rose for several months and then declined quickly. Although the early afterglow emission could not differentiate between a cocoon model or a structured jet model, the late afterglow features strongly favor the structured jet model with a large viewing angle (e.g., Lazzati et al. 2018;Troja et al. 2019), and the presence of a relativistic jet with a large viewing angle has been identified by the successful measurement of the superluminal movement of the radio afterglow emission (Mooley et al. 2018;Ghirlanda et al. 2019).
In the standard fireball model, the ejecta is assumed to be a conical outflow within which the energy density and Lorentz factor are constant (the so-called "top-hat" jet), while, beyond the ejecta, the energy density and the Lorentz factor reduce to zero abruptly (Pian et al. 2017). But in reality, the energy density and the Lorentz factor should vary with the polar angle, which has been confirmed by some numerical simulations (Aloy et al. 2005;Murguia-Berthier et al. 2017). In previous studies, structured jet models have been proposed to investigate the characteristics of the prompt and afterglow emission of LGRBs (Granot et al. 2002;Rossi et al. 2002;Zhang & Mésźaros 2002;Wei & Jin 2003), while for sGRBs it is hard to infer the energy profile of the jet since the observational data are usually rare, the exception is GRB 051221A, for which the X-ray afterglow light curves showed a flat segment and can be well accounted for by a two-component jet model (Jin et al. 2007).
Motivated by the discovery of GW170817 and GRB 170817A association, it is reasonable to assume that the sGRBs driven by binary neutron star mergers would have similar jet structure and may be observed with relative large viewing angles, so it is important to study the characteristics of these possible GW-associated sGRBs. In this paper, we investigate the observable features (such as the distribution of luminosity and viewing angles) of these off-axis sGRBs that may be associated with GW events using the structured jet model.
The paper is organized as follows. In Section 2 we describe the method we used to study the observable characteristics of nearby sGRBs, including the rate of sGRBs, the luminosity distribution and the Gaussian jet we adopted. Our results are presented in Section 3 and conclusions are given in Section 4.

The Method
When calculating the emission from a GRB jet, the most widely used jet model is a "top-hat" jet, however, the "top-hat" jet cannot explain the emission of GRB 170817A, while the structured jet model is strongly favored. For structured jet models, there is usually a critical angle θ c (the jet's core). The energy and the Lorentz factor distribution are nearly independent on the polar angle θ when θ<θ c , while for θ>θ c , the energy and the Lorentz factor drop quickly. When a GRB is viewed at an angle larger than the jet's core (θ > θ c , e.g., "offaxis"), the GRB's emission would be weak and the afterglow would rise first and then decline quickly, just as the case of GRB 170817A. The emission received by off-axis observers strongly depends on the energy density and Lorentz factor distribution within the jet. In the literature, three types of structured jet models are usually adopted: (a) the power-law distribution model q q = -  k 0 ( ) for θ>θ c and ò(θ)=ò 0 for θ<θ c (Dai & Gou 2001;Rossi et al. 2002;Zhang & Mésźaros 2002;Wei & Jin 2003) , for θ c ∼5°; and (c) the two-component jet model (Berger et al. 2003;Wu et al. 2005;Huang et al. 2006). Here we use the Gaussian distribution to describe the energy and Lorentz factor distribution (e.g., Zhang & Mésźaros 2002;Beniamini & Nakar 2019) Here, L(θ) is the luminosity per solid angle in the observer frame (Howell et al. 2019), Γ(θ) is the Lorentz factor, θ c represents the jet's core, and L 0 is the isotropic-equivalent luminosity observed on-axis. The observed gamma-ray luminosity for an observer located at angle θ obs from the jet's axis is given by (Kathirgamaraju et al. 2018;Beniamini & Nakar 2019) where δ is the Doppler factor, β(θ) is the velocity corresponding to Γ(θ), χ is the angle between the emitting material and the observer, and η γ is the efficiency of conversion from the kinetic energy to the γ-rays. It is easy to show that, for q q c obs  , the emission is dominated by "line of sight" emitters and the observed isotropic-equivalent γ-ray luminosity is just η γ ·L 0 , while, for larger viewing angles, the emission is dominated by "off line of sight" emitters and the isotropicequivalent gamma-ray luminosity would be much lower (Beniamini & Nakar 2019).
In order to study the characteristics of sGRBs, we also need to know the formation rate and luminosity function of sGRBs, which are discussed in the following subsections.

The Rate of sGRBs
It has been widely believed that sGRBs originate from the merging of binary compact objects, i.e., neutron star-neutron star (NS-NS) or neutron star-black hole (NS-BH) coalescence, and this has been confirmed by the discovery of the GW170817-GRB 170817A association. In this merging scenario, the formation rate of sGRBs can be written as (Ando 2004) where R t SF ( ) is the star formation rate (SFR), P m (t) is the probability distribution function of merging time of the binary system from its formation, t F represents the formation epoch of galaxies, and, in many papers, z(t F )=5 is usually adopted (Ando 2004). Although there remain a huge amount of uncertainties concerning the cosmic SFR history, especially in the high-redshift universe, the SFR in the low-redshift region is fairly well known by many observations at various wave bands. Among many models of SFR, two models are often used.  Figure 1 shows the SFR history for the above two models. It is obvious that the behaviors of these two models are roughly the same in the low-redshift region z<1.5, while, at z>1.5, they are significantly different. For the function SF1, it gives an exponentially decreasing formation rate, while, for the function SF2 it gives the constant formation rate at z>1.5. In this paper, we use the latest model of star formation rate (SF1) to calculate the distribution of sGRBs, and adopt the standard ∧CDM

Madau & Dickinson (Madau & Dickinson 2014)
In addition to the SFR, the merger time distribution is also essential to estimate the sGRB rate. In previous works, a simple parameterization function such as Pm(t)∝t α with a lower cutoff timescale τ is usually adopted (Porciani & Madau 2001;Schmidt et al. 2001aSchmidt et al. , 2001bGhirlanda et al. 2004), and, in general, the cutoff time t = 20 Myr m is adopted (Bulik et al. 2004;Regimbau & Hughes 2009;Meacher et al. 2015), where τ m represents the lower cutoff timescale in units of Myr. Previous studies have shown that the distribution of the merging times depends on the distribution of the initial orbital separation modeled as µ b -dN da a N , and the expected merger times follow dN/dt merge ∝τ − α , where α≡−β N /4−3/4. One can see that, even if β N varies in an extreme range from 0 to 7, α is between 0.75 and 2.5 (Ando 2004;Belczynski et al. 2017Belczynski et al. , 2018. In order to examine whether the value of α would influence the sGRB distribution significantly, we take the values of α as (−0.5, −1.0, −1.5, −2.5) to calculate the sGRB rate according to Equation (6), adopting the rate of local NS-NS mergers to be 110-3840 Gpc −3 yr −1 (Abbott et al. 2019). The results are shown in Figure 2 and we find that the different values of index α can make the sGRB rate quite different.

The Luminosity Function of sGRB
When calculating the luminosity of sGRBs, we need to know the distribution of the isotropic-equivalent luminosity L 0 . In GRB studies, the broken power-law model is most commonly used as the form of the luminosity function so, in this paper, we take the standard broken power-law distribution as follows Here L * is the critical luminosity and α L and β L are the slopes describing the low and high component of the luminosity function, respectively. Following the paper of Wanderman & Piran (2015), we take the values α L =1, β L =2, and =Ĺ 2 10 52 * erg s −1 .

k-correction
The luminosity of sGRBs can be written as    Here f(E) is the Band function of sGRBs spectra (Band et al. 1993). E min and E max denote the lower and upper cutoff values of the detector's observational energy range, α * and β * are the spectral indices, and we take the values α * =−1, β * =−2.5, and E p =800 keV (Wanderman & Piran 2015;Lien et al. 2016;Mogushi et al. 2019).

Results of Monte Carlo Simulation
Based on the above model, we can derive the distribution of observed luminosity and the viewing angle θ obs of sGRBs with the Monte Carlo simulations. By assuming the uniform distribution of the solid angle Ω(θ) in the sky, we can first calculate the Lorentz factor and the luminosity of sGRBs with Equations (1)-(5). Here we assume that L 0 satisfies the distribution of Equation (9) (Wanderman & Piran 2015) and take Γ 0 =300 (Beniamini & Nakar 2019). Several papers have investigated the emission of GRB 170817A with the structured jet model and found that θ 0 ∼5°can reproduce the observed characteristics well (e.g., Troja et al. 2019), so here we take θ 0 =5°. From Equation (6), we can obtain the redshift distribution of sGRBs. To obtain the observed luminosity and  viewing angle θ obs distribution, the selection effect must be considered. Because the Fermi-GBM telescope has a larger field of view (Burns et al. 2016), we take the Fermi-GBM 64 ms limiting flux (F lim = 2.0×10 −7 erg cm −2 s −1 ) as the detection threshold (Goldstein et al. 2017). The sGRB can be detected only when its flux is larger than this threshold.
The distribution of observed gamma-ray luminosity and viewing angle at different redshifts are shown in Figures 3 and  4. It can be seen that, for lower redshifts, the typical luminosity is relatively smaller, while the typical viewing angle is larger. For example, at z=0.045 and 0.208 (corresponding to the luminosity distance of 200 Mpc and 1 Gpc), the peak luminosities are about 8×10 47 erg s −1 and 3×10 49 erg s −1 , respectively, while the peak viewing angles are around 17°and 8°. Therefore, the luminosity/viewing-angle distribution for nearby sGRBs is quite different from previous thinking. Figure 5 gives the variation of the typical luminosity and viewing angle with redshift. It is obvious that the observed typical luminosity increases with redshift, while the typical viewing angle decreases with redshift. Especially at smaller redshifts, the luminosity rises rapidly while, at larger redshifts (such as z > 1), they increase slowly. This is because, for GRBs with larger redshifts, only an observer with smaller viewing angles (for z > 1 the viewing angle should close to the energetic core of the jet) can detect them under the structured jet scenario (see also Gupte & Bartos 2018;Howell et al. 2019). We note that our results are consistent with the distribution of the maximum viewing angles obtained by Howell et al. (2019). Figure 6 shows the luminosity and viewing-angle distributions of sGRBs within certain volumes. Since the sensitive distance of advanced LIGO/Virgo is about 200 Mpc for binary neutron star mergers or 400 Mpc for typical neutron star-black hole mergers (Abadie et al. 2010), so the sGRBs associated with present GW events are all local, at smaller redshifts. From these figures, we can see that, for sGRBs from binary neutron star mergers detected by advanced LIGO/Virgo, their luminosity mainly lie between 10 47 and 10 48 erg s −1 . For the third generation GW detector (i.e., Einstein Telescope), the redshifts of GW events can reach to ∼1. 3 In this case, the typical luminosity of sGRBs are 8×10 50 -7×10 51 erg s −1 .
We also calculate these distributions with α between −0.5 and −2.5 to see whether different models of formation rate of sGRBs will affect these distributions significantly. Figure 6 shows the variation of the typical luminosity and viewing angle with redshift for four different values of α. As can be seen from the figure, the typical luminosity of sGRBs increases very slightly with the decrease of α and, in fact, the discrepancies between these results are so small that it is hard to discriminate between different models. This is also true for the viewingangle distribution, although the typical viewing angle is somewhat smaller for smaller values of α, the difference is not significant, i.e., the distribution of viewing angles is also insensitive to the value of α. This indicates that the formation rate of sGRBs with different power-law time delay index α almost does not affect the observed luminosity and viewingangle distributions.
By comparing the total simulated number of sGRBs with what can be detected by Fermi-GBM, we obtained the detection probability of sGRBs by Fermi-GBM; the results are shown in the left panel of Figure 7. We can see that the detection probability decreases with redshift; when the redshift is between 0 and 0.2, the detection probability drops very sharply and, at z=1, the detection probability is only about 10 −3 . We also give the normalized distribution of number of observed sGRBs at different redshifts, as shown in the right panel of Figure 7. Here we have calculated the distribution of number of sGRBs for four different values of α and it is very interesting that these distributions are obviously different. For smaller values of α, the redshift corresponding to the peak number of observed sGRBs is larger, and at redshift z=1, the observed number of sGRBs for α=−2.5 is much larger than that for α=−0.5, −1, and −1.5, thus it provides a better way to determine the values of α.
We use Bayesian method to estimate the number of sGRBs required to distinguish different α values. The posterior of α after observing a series of events with redshift can be obtained under the standard Bayesian framework: a a a µ P D L D P 12 In Equation (12), a L D ( |ˆ) is the likelihood of the data series D predicted by the theoretical distribution calculated from a specific α ( a N z *( ), as shown in Figure 7). a P (ˆ) is the prior, which is taken to be uniform in this study, and the likelihood term can be given by Equation (13). We draw random redshift samples with z<1 from a N z *( ) for different α values (−0.5, −1, −1.5, and −2.5), then recover the α values by calculating their posterior distributions described by Equation (12). The number of samples progressively increases, so that we can investigate how the constraint on α changes with the observed GRB numbers. This process is repeated 100 times for each injected α, and the averaged 95% confidence intervals of the recovered α are shown in Figure 8. The number of samples needed to let the 95% confidence intervals lie between the two adjacent α with respect to the injected one are 220, 290, 345, and 540 for a=−0.5, −1.0, −1.5, and −2.5, respectively.
where R sGRB (z) is the formation rate of sGRBs as given by Equation (6) and f sGRB (z) is the detection probability as shown in Figure 7. We use the total time-averaged observable sky fraction of Fermi-GBM which is given as Λ GRB =0.60 (Burns et al. 2016). In the top-hat model, the conversion of the observed rate to the intrinsic rate typically needs a geometrical beaming factor f b , and where θ j is the jet half-opening angle. However, for the structure jet model, the effect of beaming factor has actually been included in the detection probability. Then we use Equation (14) Table 1. As can be seen from them, the change of a has little effect on the observed number of sGRBs at low redshift; this is not surprising since, from Figure 2, we can see that the changing a has negligible effect on the formation rate of sGRBs at low redshift.  (Howell et al. 2019). By comparing this with our results, we find that the smaller value of a (such as a=−2.5) is disfavored.
We also compare our results of gamma-ray luminosity distribution with sGRBs with known redshifts detected by Fermi-GBM, as shown in Figure 10 (the data is in Table 2). It is evident that our results are consistent with the observation. In Figure 11, We perform the Kolmogorov-Smirnov test for the   Figure 10. The relation between gamma-ray luminosity and redshift. The red circles represent our simulation results, the green triangle represents sGRB 170817A, and the blue circles are the sGRBs with known redshifts observed by Fermi-GBM (in Table 2). The error bars are 1σ errors. cumulative luminosity distribution of our simulated sGRBs and those with known redshifts. The chance probabilities is 0.48, which indicates that our simulation results can reproduce the observed luminosity distribution well.

Discussion and Conclusion
The joint detection of GW and gamma-ray emission from the binary neutron star merger event GW170817/GRB 170817A marks a breakthrough in the field of multi-messenger astronomy and confirms the association of sGRBs with binary neutron star mergers. The peculiar behavior of the prompt and afterglow emission of GRB 170817A indicates that this nearby sGRB was seen from a structured jet with large viewing angle, and this speculation has been identified by the successful measurement of the superluminal movement of the radio afterglow emission (Abbott et al. 2017a(Abbott et al. , 2017bGoldstein et al. 2017;Savchenko et al. 2017).
Since the beginning of the third observation run of advanced LIGO/Virgo, a number of NS-NS and NS-BH merger events/ candidates have been detected, such as the notable S190425z (LIGO Scientific Collaboration & Virgo Collaboration 2019a, 2019b) and S190426c (LIGO Scientific Collaboration & Virgo Collaboration 2019c). Although no sGRBs associated with them have been detected, it is consistent with our prediction about the expected detection rate of nearby sGRBs (see Figure 9). The detection of these NS-NS and NS-BH merger events imply that a number of potential sGRBs really exist within the sensitive distance of advanced LGO/Virgo and they could probably be detected by the next generation of more sensitive gamma-ray detectors, so it is important to investigate the properties of these nearby sGRBs. In this paper, we use Monte Carlo simulation to calculate the gamma-ray luminosity and viewing-angle distribution of nearby sGRBs based on the structured jet model, and obtain the evolution of these distributions with redshift. We find that, for sGRBs from binary neutron star mergers detected by advanced LIGO/Virgo, their luminosity mainly lies between 10 47 and 10 48 erg s −1 .
In our calculation, we assume the formation rate of sGRBs is proportional to the SFR. Although the SFR is uncertain at high redshifts, it is well known at low redshifts, and various SFR models are only slightly different at low-redshift region z<1.5, so for the nearby sGRBs that we are most interested in, the uncertainty of SFR have nearly no effects on our results. Here we use the latest SFR model of Madau & Dickinson (2014) to calculate the sGRB rate. In addition to the SFR, the merger time distribution is also essential to estimate the sGRB rate. From Figure 6, we see that the typical luminosity and viewing angles of sGRBs seem to be insensitive to the value of α, which means that we cannot identify the value of α through observed luminosity and viewing-angle distribution. However, it is very interesting that, for different values of α, the distribution of numbers of observed sGRBs are quite different, so it is possible to determine the value of α through observed distribution of the number of sGRBs. We used the Bayesian method to make a quantitative analysis and found that the value of α may be identified when the number of observed sGRBs with known redshifts is more than 200.
During the advanced LIGO/Virgo O3 run a number of NS-NS and NS-BH merger events/candidates have been detected, but no corresponding sGRBs have been observed. This is consistent with our results that the expected detection rate of sGRBs is only about 1 yr −1 by Fermi-GBM within the distance of several hundred Mpc, so developing more sensitive GRB detectors is very important for future studies of GRBs, especially for nearby sGRBs that may be weak and associated with GW events. The GECAM (Gravitational Wave Electromagnetic Counterpart All-sky Monitor) mission is the planned Chinese space telescope to monitor the GRBs coincident with GW events, which will launch in late 2020. It has a FOV of 100% all-sky and it has lower detection threshold than Fermi-GBM, about the same as Swift (Zhang et al. 2019;Zheng 2019). Our calculation shows that the peak luminosity is about 6×10 46 erg s −1 and 4.0×10 48 erg s −1 at luminosity distance of 200 Mpc and 1 Gpc, respectively, and the number of sGRBs that can be detected by GECAM are 1.04−35.96 yr −1 and 16.97−593.91 yr −1 within 200 Mpc and 1Gpc, respectively, which is about five times larger than that can be detected by Fermi-GBM.