GRBs luminosity function synthesized from Swift/BAT, Fermi/GBM and Konus-Wind data

Gamma-ray bursts (GRBs) are extremely powerful explosions that involve either the death of massive stars or the merger of compact stellar objects. The equivalent isotropic energy of GRBs can exceed $10^{54}$ erg, and their nonthermal spectra peak between 10 and $10^{4}$ keV (e.g. Atteia et al. (2017)). Their light curves display intense pulses that can last anywhere from less than a second to hundreds of seconds, they are often irregular and sometimes include multiple pulses. Although the radiation emitted by GRBs is believed to come out mainly through jets (with collisions inside), the precise formation mechanism behind these jets is not fully understood (e.g. Le and Mehta (2017)).

Given the complex nature of their occurrence and unfolding, it is not surprising that GRBs are among the most studied topics and phenomena in astrophysics today, both theoretically and observationally. Many theoretical models have been proposed to try and explain this very spectacular phenomenon; but, unsurprisingly given the astounding amount of energy involved and the short time scales of these explosions, none of these models has been fully successful ((Dainotti et al., 2019)).

There are currently many observation missions, both from the ground and from space, that are involved in the study of GRBs in the different bands of the electromagnetic spectrum. Furthermore, a new observational tool has recently emerged, namely gravitational waves, which are produced by the merger of neutron stars or black holes. This new tool is already providing insight into the formation mechanisms of short, merger-type GRBs and has become an integral part of the multi-messenger approach to studying GRBs (Abbott and al. 2017, Drago et al. 2018, Mészáros 2019, and others).

The study of the luminosity function (LF), which is defined as the number of bursts per unit interval of luminosity, is part of the verification and validation of the theoretical models of the physical processes involved in the occurrence of GRBs. It allows us, among other things, to predict the production rates of GRBs and the probability of detecting gravitational waves associated with some of them. Furthermore, it can provide insights into the star-formation rate.

Previous research on luminosity functions has been done on different data samples, most of which contained a limited number of gamma-ray bursts (Schmidt, 1999, 2001; Sethi and Bhargavi, 2001; Yonetoku et al., 2004; Liang et al., 2007; Virgili et al., 2009; Butler et al., 2010; Cao et al., 2011; Salvaterra et al., 2012; Wanderman and Piran, 2015; Amaral-Rogers et al., 2016; Paul, 2017; Kinugawa et al., 2019; Lan et al., 2019). Our present work uses three larger data samples from three satellites/instruments, namely Swift/BAT (251 GRBs), Konus/Wind (152 GRBs), and Fermi/BAT (37 GRBs), and then enlarges the data set with simulated (artificial) GRBs.

In Section 2, we describe how we select our samples. Sections 3 and 4 provide details of how we calculate the luminosities and the luminosity functions (for each sample), respectively. In Section 5, We present our results and then discuss them in Section 6, ending with concluding remarks and a mention of future research prospects in Section 7.

For the Swift/BAT and Fermi/GBM samples, the luminosity at the peak of the flux can be calculated using the spectral data of each burst, characterized by a known redshift. The 1-second photon flux at the peak can be found in the Swift/BAT and data Fermi/GBM. This flux gives the total number of photons per unit area and per unit time, regardless of their individual energies. The Swift/BAT energy spectrum is divided into four bands: 15-25 keV, 25-50 keV, 50-100 keV, 100-350 keV. The Fermi/GBM covers the wider band 10-1000 keV. In this work, we use a cut-off power-law (Band et al., 1992), which is characterized by two spectral parameters: the observed peak energy, $E^{obs}_{p}$, and the spectral index alpha ($\alpha$):

$$N_{E}(E)=AE^{\alpha}~{}e^{-E/E_{0}},$$

(1)

where $E_{0}=E^{obs}_{p}/(2+\alpha)$.

The bolometric luminosity of the 1-second isotropic peak flux, denoted by $L_{p}$, is the maximum energy radiated per unit time in all space. It is calculated by integrating the $EN_{E}$ function in the energy band corresponding to the observed gamma radiation band in the source’s frame, i.e. $E_{1}$ = 1 keV to $E_{2}$ = $10^{4}$ keV. And because of cosmological effects, the corresponding observed energy band is: $E_{1}/(1+z)$ to $E_{2}/(1+z)$.
Thus, the k-corrected $L$ is calculated via:

$$L=4~{}\pi~{}d_{L}^{2}~{}P_{\gamma}~{}k_{c}.$$

(2)

Here $L$ is k-corrected with the method developed by (Bloom et al., 2001), and $P_{\gamma}$ is the peak flux expressed as $ph~{}cm^{-2}~{}s^{-1}$. The k-correction factor is defined by (Yonetoku et al., 2004; Rossi et al., 2008; Elliott et al., 2012) as:

$$k_{c}=\frac{\int_{E_{1}/(1+z)}^{E_{2}/(1+z)}EN(E)dE}{\int_{E_{min}}^{E_{max}}N(E)dE}\ .$$

(3)

The cosmological distance $d_{L}$ is given by the equation:

$$d_{L}=\frac{(1+z)c}{H_{0}}\int_{0}^{z}\frac{dz^{\prime}}{\sqrt{\Omega_{M}(1+z^{\prime})^{3}+\Omega_{\Lambda}}}\ ,$$

(4)

where we adopt the following cosmological parameters: $\Omega_{M}$ = 0.27, $\Omega_{\Lambda}$ = 0.73, and $H_{0}$ = 70 $\mathrm{km/s/Mpc}$ (Komatsu et al., 2009).

These integrals are performed numerically using the time-resolved spectral parameters given by Swift/BAT and Fermi/GBM.

With this data, as explained above, we calculate the luminosity of each burst, taking into account the error on each quantity used. The results of this calculation are presented in Tables LABEL:tabSwift and LABEL:tabFermi. The errors are evaluated using the Monte Carlo method. The quantities needed in the calculation of the luminosity ($E_{p}$, $\alpha$ , flux) are assumed to follow a normal distribution. For each GRB, by making N random draws from these distributions, we obtain the two extremum values of the luminosity. The subtraction of these two values from the calculated value without uncertainty gives the errors.

For the Konus Wind sample, given the difficulty of obtaining the spectral parameters ($\alpha$, $P_{\gamma}$) for all the bursts, we preferred to use the correlation relation of Yonetoku et al. (2004), which relates the luminosity to $E_{p}$. We used the data for z and $E_{p}$ given in Minaev and Pozanenko (2019). We also used the improved Yonetoku et al. (2004) correlation relation:

$$\frac{L}{10^{52}}=(2.34^{+2.29}_{-1.76})\times 10^{-5}\bigg{[}\frac{E_{p}(1+z)}{1~{}keV}\bigg{]}^{2.0\pm 0.2}\ .$$

(5)

In our Monte Carlo simulation, we draw values of $L$ using the following relation (which is the same as (5)):

$$\frac{L(j)}{10^{52}}=\alpha\times 10^{-5}\bigg{[}\frac{E_{p,j}(1+z)}{1~{}keV}\bigg{]}^{q},$$

(6)

where $L(j)$ and $E_{p,j}$ are the luminosity and the energy corresponding to the jth GRB in our data table. To estimate the error on $L(j)$, we assume that $\alpha$ and $q$ are random quantities obeying a normal distribution, such that: $\alpha=2.34^{+2.29}_{-1.75}$ and $q=2.0\pm 0.2$. We also assume that the peak energy $E_{p}(j)$ is normally distributed. The results of this calculation are presented in Table LABEL:tabKW.

The total number of gamma-ray bursts observed by a telescope within an aperture solid angle $\Omega$ and during an effective period of of activity T (in years) is:

$$N_{T}=\frac{\Omega T}{4\pi}\int_{0}^{z}\frac{R_{GRB}(z)}{1+z}\frac{dV}{dz}dz\int_{L_{min}}^{L_{max}}\varphi(L)dL,$$

(7)

where $\varphi(L)$ is the luminosity function and $R_{GRB}(z)$ is the rate of occurrence of GRBs at redshift z expressed in $Mpc^{-3}yr^{-1}$. The factor (1+z) introduces the cosmological expansion factor and $\frac{dV}{dz}$ is the variation of the volume of the universe as a function of the redshift z, given by the equation :

$$\frac{dV}{dz}=\frac{c}{H_{0}}\frac{4\pi d_{L}^{2}}{(1+z)^{2}\sqrt{\Omega_{M}(1+z)^{3}+\Omega_{\Lambda}}}.$$

(8)

The GRB luminosity function is usually defined as the number of bursts per unit of luminosity (Sethi and Bhargavi, 2001; Daigne et al., 2006; Zitouni et al., 2008; Cao et al., 2011; Amaral-Rogers et al., 2016; Paul, 2017):

$$\varphi(L)=\frac{dN}{dL}.$$

(9)

To predict the total number $N_{T}$ of gamma-ray bursts that will be observed by a telescope, we need to know the rate of occurrence of GRBs and their luminosity function $\varphi(L)$. The $R_{GRB}$ rate is very poorly known despite the fact that there are several attempts based on broad assumptions.

In the present work, we assume that the rate of long gamma-ray bursts is proportional to the star formation rate (SFR). Indeed, it is now accepted that long gamma-ray bursts are due to the collapse of massive stars, thus the burst rate should be proportional to the supernova rate through a constant k parameter (Porciani and Madau, 2001; Daigne et al., 2006):

$$R_{GRB}=k\times R_{SN},$$

(10)

where $R_{SN}$ is the Type II supernova comoving rate. Moreover, these authors assume that the rate of supernovae is proportional to the rate of massive star formation:

$$R_{SN}=0.0122\ M_{\odot}^{-1}\times R_{SFR},$$

(11)

where $R_{SFR}$ is the star formation rate, which is determined by observations only up to a redshift $z\simeq 2-2.5$. Beyond this value, several hypotheses are made in the literature, from which we retain the three most common ones: $SFR_{1}$, where the star formation rate decreases beyond $z\geq 2$; $SFR_{2}$, where the rate remains constant beyond $z\geq 2$; and $SFR_{3}$, where the rate continues to increase beyond $z\geq 2$ (Daigne et al., 2006). These three hypotheses are expressed by the following relations (Porciani and Madau, 2001):

$$SFR_{1}(z)=0.323\ h_{70}\ \frac{e^{3.4z}}{e^{3.8z}+45}\quad[M_{\odot}~{}yr^{-1}Mpc^{-3}]$$

(12)

$$SFR_{2}(z)=0.160\ h_{70}\ \frac{e^{3.4z}}{e^{3.4z}+22}\quad[M_{\odot}~{}yr^{-1}Mpc^{-3}]$$

(13)

$$SFR_{3}(z)=0.215\ h_{70}\ \frac{e^{3.05z-0.4}}{e^{2.93z}+15}\quad[M_{\odot}~{}yr^{-1}Mpc^{-3}]$$

(14)

In our work, we have found it useful to use an equivalent form of the luminosity function presented by Butler et al. (2010). This form can be considered as a scale-free function (Sethi and Bhargavi, 2001). Moreover, we are interested in determining the parameters (detailed below) that characterize this luminosity function, so we will use a form of normalized luminosity function defined as follows:

$$\phi(L)=\frac{1}{N_{T}\times L}\frac{dN}{d\log{L}}.$$

(15)

To calculate the GRB luminosity function, we use the samples of GRBs (and their data) from the Swift/BAT, the Fermi/GBM, and the Konus/Wind satellites/instruments. With the exception of the redshift, all the parameters used in our calculations are given with their respective uncertainties. We have assumed that each parameter value is a normally distributed quantity characterized by the density function:

$$f(x)=\frac{1}{\sigma\sqrt{2\pi}}\exp{-\frac{(x-\mu)^{2}}{2\sigma^{2}}},$$

(16)

where $\mu$ represents the mean value of the measured quantity, and $\sigma$ corresponds to its uncertainty, both given in Table 8. For example, for each GRB, the energy at the peak $E_{p}$ given in the table corresponds to $\mu$ and the $\delta E_{p}$ corresponds to $\sigma$.

Let us now compare and contrast our work and results to those of other researchers.

Like our own, the investigation by Wanderman and Piran (2015) involved a joint analysis of GRBs from different satellites, with data from Swift, Fermi, and CGRO-BATSE, aiming to determine the luminosity function and the formation rate of GRBs. However, unlike our study, their investigation was restricted to short GRBs, and more specifically to non-collapsar short bursts, having found that their short-GRBs sample was ”contaminated” by collapsars, which affected their data sample. With those restrictions, they obtained a luminosity function that was best fit with a broken power law with a break at $L_{b}=2\times 10^{52}erg.s^{-1}$, which is consistent with our results.

The investigation by Pescalli, A. et al. (2016) is most relevant to us because, like our study, theirs examines the luminosity function of GRBs and briefly explores whether GRBs undergo luminosity evolution. However, these authors restrict their data sample to long Swift GRBs, which reduces their sample to only 99 bursts compared to the 439 bursts in ours. Furthermore, our data is augmented by thousands of simulated ”artificial” bursts for simulation purposes. Despite these differences, Pescalli, A. et al. (2016) find that their data can be best fit with a broken power law with $\log{(L_{b})}=51.45$, which is consistent with what we obtain. However, although their high luminosity index $\nu_{2}(=-1.84)$ is consistent with our value, their low luminosity index $\nu_{1}(=-1.32)$ is not, indicating that the low-end of their luminosity function is steeper than ours. The difference could be due to their data sample being restricted to long GRBs taken from Swift only, whereas our sample consists of data combined three sub-samples. We should note that a close examination of our results reveals that among our three sub-samples, the $\nu_{1}$ value for the Swift sub-sample was significantly the steepest among the three, which is at least qualitatively in agreement with Pescalli, A. et al. (2016), and it might be hinting that the lower-end of the luminosity function is steeper for the Swift GRBs in particular. Moreover, like us, Pescalli, A. et al. (2016) found that there was a correlation between the luminosity and the redshift but concluded that this evolution should not be interpreted as being purely due to luminosity evolution. We agree with them, as the correlation with the redshift could be due to density evolution as well, and like Pescalli, A. et al. (2016) we admit that this issue needs further work, beyond the scope of our current study, because simulating both the luminosity and density evolutions and extricating the impact of each on the luminosity function and finding the extent and shape of its evolution is a separate and full-fledged topic. This is also the conclusion reached by Lan et al. (2019) who analyzed the same data sample used by Pescalli, A. et al. (2016) to investigate the GRB luminosity function and utilize it to determine the star-formation rate.

The investigation by Amaral-Rogers et al. (2016) takes a distinct approach from ours. These authors explore the pulse luminosity function for Swift GRBs. They use a sample of 118 bursts to generate a total of 607 pulses and then compute the pulse luminosity function by employing three different GRB formation rate models. Despite this markedly different approach that includes contributions from all fitted pulses and not just from the brightest ones, the best-fit parameters for their luminosity function, $L_{b}=1.5\times 10^{51}erg.s^{-1}$, $\nu_{1}=-0.70$, and $\nu_{2}=-1.69$, are consistent with our results, although once again their low-end luminosity function is somewhat steeper than what we obtain.

The investigation of the GRB luminosity function is an important issue in GRB research because it gives insight into the physics involved in the production mechanism of GRBs. Not only does it shed light on the range of luminosities that GRBs can have, it also reveals whether there is a ”scale” for GRB luminosities. The fact that all three sub-samples that we explored had a broken power-law luminosity function, with break luminosities consistent with one another, strongly indicates that the physics behind the formation of GRBs is not ”scaleless”, a result which can be useful when probing other physical aspects of GRBs.

In this study, we performed separate analyses of Swift/BAT, Fermi/GBM, and Konus/Wind data for all GRBs with known redshifts, in the aim of determining their luminosity functions. In doing that, we left out the formation rate of the GRBs from our considerations. We obtained the best parameters for the luminosity functions using the likelihood method. We found that they can be well represented by a broken power law -the so-called broken luminosity function.

The luminosity function expressed in log-log is well fit with a broken power-law with a break $L_{b}$ in the luminosity, which varies from $1.2\times 10^{52}~{}erg.s^{-1}$ to $6.3\times 10^{52}erg.s^{-1}$. The exponent-index of the low luminosity part, $\nu_{1}$, varies between -0.23 and -0.47 and that of the high luminosity part, $\nu_{2}$, varies between -1.79 and -2.68. These power law indices are consistent with those found in previous studies (values and references are given in Table 9).

In the future, we plan to look at both the density and luminosity evolution of GRBs, which is something that we only brushed briefly in this study. Such an investigation is important because it will cast light on several critical issues such as the star-formation rate, the degree to which GRBs can be utilized as ”standard candle”, and the evolution of the luminosity function itself.

The authors gratefully acknowledge the use of the online Swift/BAT table compiled by Taka Sakamoto and Scott D. Barthelmy and Evans, P. A. et al. (2010). The authors also aknowledge the use of the online Fermi/GMB table compiled by von Kienlin et al. (2014), Gruber et al. (2014) and Bhat et al. (2016). We thank the General Directorate of Scientific Research and Technological Development (DGRSDT), Algiers, Algeria for the financial support.