Abstract

Gamma-ray bursts (GRBs) are extremely powerful explosions that have been traditionally classified into two categories: long bursts (LGRBs) with an observed duration T₉₀ > 2 s, and short bursts (SGRBs) with an observed duration T₉₀ < 2 s, where T₉₀ is the time interval during which 90% of the fluence is detected. LGRBs are believed to emanate from the core-collapse of massive stars, while SGRBs are believed to result from the merging of two compact objects, like two neutron stars. Because LGRBs are produced by the violent death of massive stars, we expect that their redshift distribution should trace the star-formation rate (SFR). The purpose of our study is to investigate the extent to which the redshift distribution of LGRBs follows and reflects the SFR. We use a sample of 370 LGRBs taken from the Swift catalog, and we investigate different models for the LGRB redshift distribution. We also carry out Monte Carlo simulations to check the consistency of our results. Our results indicate that the SFR can describe the LGRB redshift distribution well for high redshift bursts, but it needs an evolution term to fit the distribution well at low redshift.

Keywords

Gamma-Ray Bursts, Redshift Distribution, Star-Formation Rate

Share and Cite:

1. Introduction

Although more than 50 years have passed since the first paper on gamma-ray bursts (GRBs) was published [1] , we still have many questions regarding their nature and characteristics. They are the focus of many research studies due to their potential as cosmological probes [2] - [8] . GRBs are traditionally classified into two main types: long GRBs (LGRBs) with T₉₀ > 2 s and short GRBs (SGRBs) with T₉₀ < 2 s [9] , where T₉₀ is the observed duration during which 90% of the fluence is detected. It is theorized that the origin of LGRBs is the core-collapse of massive stars, while SGRBs are thought to originate from compact object mergers [10] [11] [12] . Hence, it is expected that there will be a correlation between the star formation rate (SFR) and the LGRB formation rate, which encouraged researchers to investigate the redshift distribution of GRBs, and especially that of LGRBs.

However, it is important to note that the connection between the SFR and the LGRB redshift distribution is a controversial one, with many conflicting results in the literature. The controversy lies in whether there is an observed excess of GRBs at low or high redshift compared to what is implied by the SFR. The study by [13] is one of the first major contributions to this problem. Although the connection between the SFR and LGRB redshift distribution had been studied before, see for example [14] , the study by [13] is one of the first to account for completeness and to place constraints on the GRB luminosity function and its evolution. Their study found evidence for a possible luminosity or density evolution of LGRBs compared to what is expected from the SFR [13] . A recent study, but with a larger data sample, by [15] found similar results. An earlier study by [16] , which utilized 127 LGRBs observed by the Swift satellite, found an excess in GRBs at low redshift compared to the SFR. Although some studies backed the low-redshift excess claims [17] [18] [19] [20] , other studies found an excess in GRBs, not at low redshift, but at high redshift [21] [22] .

Moreover, many studies have opposed the excess GRB accounts, especially at low redshifts. The study by [23] was one of the earliest to do so, where they carried out Monte Carlo simulations and compared their results with previously published claims of low redshift excess – the aforementioned [16] [19] studies. They showed that the excess GRBs result from the incompleteness of the dataset, and that a complete dataset shows no excess in GRBs. The study by [23] also conducted a non-parametric C⁻ investigation. However, their non-parametric study seems misleading as discussed by [24] who made an in-depth study about the non-parametric methods and how their misuse can lead to incorrect results. They analyzed datasets of previous GRB luminosity/redshift evolution studies, like the previously mentioned [16] [18] [20] [23] . They found that many studies underestimated the detection threshold of GRB satellites which resulted in incomplete datasets, and hence showed an excess in GRBs. They also created a Monte Carlo simulation, which they found to agree with the Monte Carlo simulations conducted by [23] ; in other words, excess GRBs resulted from incomplete datasets. The investigation by [25] also found no excess GRBs at any redshift.

More recently, the excess GRBs at low redshift claims have started to reemerge. They are mainly driven by studies conducted by [26] [27] . In the first study, [26] formed two sets of GRBs composed of 54 and 79 LGRBs, with 61% and 67% completeness, respectively. They found an excess of LGRBs at low redshift, even with the completeness of the dataset taken into consideration. The second study [27] focused on potential reasons for this excess of LGRBs. They split their datasets into low and high luminosity LGRBs. They found that the high luminosity LGRB dataset follows the SFR without any excess, while the low luminosity LGRB dataset did show an excess of LGRBs at low redshifts. They argue that their results may hint at new origins for low luminosity LGRBs, thus disconnecting the LGRB formation rate from the SFR. Interestingly, there were some recent reports that suggest that LGRBs can be produced by compact object mergers, instead of massive star collapse [28] [29] . Moreover, there have been studies that investigated the effect of metallicity on the LGRB distribution and have found it to be influential [30] [31] . Perhaps these can be potential reasons for the excess GRB dilemma. Although as it stands, there is no clear answer to this problem.

The aim of our current study is to investigate the distribution of LGRBs with redshift and to fit it with several proposed models. A full LGRB data set is used without constraints, including biases in detection, to study how it is linked to the star formation rate. We aim to examine the nature of the evolution term and see how it impacts the redshift distribution of LGRBs. Section 2 provides a detailed description of our method, while Section 3 presents our results and analysis. We provide a brief conclusion in the last section.

2. Method

The data sample that we use is collected from the Swift catalog¹. The data sample is composed of 370 LGRBs with T₉₀ > 2 s. The dataset is found in Table A1 in Appendix A. We divided the sample into 30 bins for $0\le z\le 10$ with equal bin

sizes $dz=\frac{1}{3}$ , and we then normalized the distribution. Mathematically, the distribution of GRBs can be written as:

$\Phi \left(z\right)=\frac{\text{d}N}{\text{d}z}=\frac{\text{d}N}{\text{d}V}\frac{\text{d}V}{\text{d}z}=A\frac{{\psi }_{*}\left(z\right)}{1+z}\varphi \left(z\right)\frac{\text{d}V\left(z\right)}{\text{d}z}$ (1)

where ${\psi }_{*}\left(z\right)$ is the star-formation rate density (SFRD), $\varphi \left(z\right)$ is a term that contains all contributions to the distribution other than the SFRD, $\frac{\text{d}V}{\text{d}z}$ is the

co-moving volume element and A is a normalization constant. In this paper, we do not assign any specific “meaning” to $\varphi \left(z\right)$ and A. All possible contributions such as the luminosity function, density evolution, sensitivity of the detectors and others are included in $\varphi \left(z\right)$ . Our aim is to gain insight on the nature of $\varphi \left(z\right)$ without presuming any prior connections or contributions from other sources. This is done to focus our work on the link to the SFRD, and to evaluate how significant the other contributions are to the redshift distribution of LGRBs. We caution here that since detection bias is also considered in Equation (1), then the equation represents the GRB detection distribution, rather than the formation distribution.

The SFRD that is typically used in the literature in units of ${\text{M}}_{\odot }\cdot {\text{yr}}^{-1}\cdot {\text{Mpc}}^{-3}$ is [32] [33] :

${\psi }_{*}\left(z\right)=\frac{0.0157+0.118z}{1+{\left(\frac{z}{3.23}\right)}^{4.66}}$ (2)

We adopt the ΛCMD model in this study, which gives:

$\frac{\text{d}V}{\text{d}z}=\frac{4\pi c{D}_L^2}{H_0{\left(1+z\right)}^2\sqrt{{\Omega }_{\text{M}}{\left(1+z\right)}^3+{\Omega }_{\text{Λ}}}}$ (3)

where $D_L$ is the luminosity distance, which is given by:

$D_L\left(z\right)=\frac{\left(1+z\right)c}{H_0}{\displaystyle {\int }_0^z\frac{\text{d}z}{\sqrt{{\Omega }_{\text{M}}{\left(1+z\right)}^3+{\Omega }_{\text{Λ}}}}}$ (4)

and $H_0=70.8\text{km}/\text{s}/\text{MPc}$ , ${\Omega }_{\text{M}}=0.27$ , and ${\Omega }_{\text{Λ}}=1-{\Omega }_{\text{M}}=0.73$ are the adopted cosmological parameters.

We investigated and fitted our data sample with multiple models and the goodness of the fit for each model was determined using the least ${\chi }^2$ test:

${\chi }^2=\underset{i}{{\displaystyle \sum }}\frac{{\left(N_i^{\text{swift}}-N_i^{\text{cal}}\right)}^2}{{\sigma }_i^2}$ (5)

where ${\sigma }_i=\sqrt{N_i^{\text{swift}}}$ is the Poisson error (68% confidence). The best fitting model is the one that gave the least ${\chi }^2$ value.

The first model that we tried was a “no contributions” model, $\varphi \left(z\right)=\text{constant}$ , with the SFRD and the cosmological model being the only contributors. Then, we tested five different contribution terms $\varphi \left(z\right)$ : a power law, a broken power law, a triple power law, an exponential term, and an exponential-power law, as follows:

${\varphi }_{\text{pow}}\left(z\right)={\left(1+z\right)}^{\alpha }$ (6)

${\varphi }_{\text{broken}}\left(z\right)=\{\begin{array}{c}\frac{{\left(1+z\right)}^{\alpha }}{{\left(1+r_1\right)}^{\alpha }}\text{for}z\le r_1\\ \frac{{\left(1+z\right)}^{\beta }}{{\left(1+r_1\right)}^{\beta }}\text{for}z>r_1\end{array}$ (7)

${\varphi }_{\text{triple}}\left(z\right)=\{\begin{array}{ll}\frac{{\left(1+z\right)}^{\alpha }}{{\left(1+r_1\right)}^{\alpha }}\hfill & \text{for}z\le r_1\hfill \\ \frac{{\left(1+z\right)}^{\beta }}{{\left(1+r_1\right)}^{\beta }}\hfill & \text{for}{r}_1<z\le r_2\hfill \\ \frac{{\left(1+r_2\right)}^{\beta }}{{\left(1+r_1\right)}^{{\delta }_2}}\frac{{\left(1+z\right)}^{\gamma }}{{\left(1+r_2\right)}^{\gamma }}\hfill & \text{for}{r}_2\le z\hfill \end{array}$ (8)

${\varphi }_{\exp }\left(z\right)={\text{e}}^{\mu z}$ (9)

${\varphi }_{\exp \text{-pow}}\left(z\right)={\left(1+z\right)}^{\alpha }{\text{e}}^{\mu z}$ (10)

where $\alpha ,\beta ,\gamma ,\mu ,r_1$ , and $r_2$ are the optimization parameters². Although we do not force any specific form on $\varphi \left(z\right)$ , the expressions that we selected and chose to investigate are based on our intuitive expectations of the insight that they might provide. The broken and triple power laws were selected to give insight on how the distribution behaves at different redshift ranges. The power and exponential laws were selected to be compared with the broken/triple power law to see if the evolution is region specific, or generic to the whole distribution.

In this work, we chose two optimization methods. The first was a fitting method. The fitting was done for $\alpha ,\beta ,\gamma ,\mu$ , while $r_1$ and $r_2$ were found by looping over them and finding the optimal values that minimize ${\chi }^2$ . The second method involved using the Metropolis-Hastings Markov Chain Monte Carlo (MCMC) algorithm [34] . The method used is described by Foreman-Mackey and Hogg and utilizes their Python libraries emcee³ and corner.py⁴ [35] [36] . The log of the likelihood function that we wish to maximize is of the form:

$\ln \mathcal{L}=-\frac{1}{2}\underset{i}{{\displaystyle \sum }}\left[\frac{{\left(N_i^{\text{swift}}-N_i^{\text{cal}}\right)}^2}{{\sigma }_i^2}\right]=-\frac{1}{2}{\chi }^2.$ (11)

Then ${\chi }^2$ can be calculated as follows:

${\chi }^2=-2\ln \mathcal{L}$ (12)

Here maximizing the likelihood function is equivalent to minimizing the ${\chi }^2$ value. However, since the models used include different numbers of parameters, the Akaike information criterion (AIC) was used to assess the goodness of the fits and to compare the different models [37] [38] :

$\text{AIC}=2k+{\chi }^2=2k-2\ln \left(\mathcal{L}\right)$ (13)

where k is the number of parameters. The better fitting model is the one that minimizes the AIC value.

We then used the MCMC simulations to fit the data to the following:

1) A new SFRD-like function—call it the GRB detection density function (GRB-DD):

${\psi }_{\text{GRB}}\left(z\right)=\frac{\alpha +\beta z}{1+{\left(\frac{z}{\gamma }\right)}^{\delta }}$ (14)

2) A broken power law (see Equation (7)).

3) A triple power law (see Equation (8)).

4) An exponential-power law (see Equation (10)).

For each case, we report the mean of the probability distribution and its deviation, the best fitting value, the ${\chi }^2$ value, and the AIC.

3. Results

The normalized redshift distribution of the data sample used is shown in Figure 1(a).

The main highlight of the results is the immediate spike of the number density, peaking at $z~1.25$ . Comparing that with the density calculated with no evolution (Figure 1(b)), we can see that the “no evolution” density peaks at $z=2$ with a subtle increase compared to the Swift data. Focusing on the low and high redshift regions, we notice the calculated density matches the Swift distribution at high redshifts. However, at low redshifts it is a shifted version of the Swift data, shifted to the right by about $\text{Δ}z~0.5$ . To solve this issue, one might shift the density left but that would not solve the problem, as then it would underestimate the number density at $z\gtrsim 4$ , which presents the main issue with this model. We need a term $\varphi \left(z\right)$ to be included in the density that scales it at low redshift while maintaining a good fit at high redshifts, like the SFRD (Figure 1(b)).

3.1. Fitting Results

We fitted the data sample with the proposed functions (Equations (6) to (10)), and we extracted the relevant parameters. The results are shown in Table 1. Based on the AIC test, the best fit occurred when we used the triple power law, the second best fit when we used the broken power law, and the third best fit when we used the exponential-power law. The power law and the exponential cases show very poor fitting. These results highlight the need for an additional $\varphi \left(z\right)$ factor and suggest that it is redshift “region” specific.

3.2. MCMC Results

The MCMC results are shown in Table 2. It is notable that the results agree mostly with the fitting results found (within 1σ) and lead to the same conclusion. Thus, we can conclude that the results are consistent with each other. All proposed models show good fitting, with the triple power law showing the best fit. Figure 2 plots all the models for comparison.

Now focusing on the variation between the fitting of the contributions at low and high redshifts, we notice that the results show a significant difference between the two redshift regions. This is clearly exhibited by the fitting and MCMC results, particularly in the broken power law and the triple power law cases. These cases mainly split the distribution into two regions, separated at around $z~1.75$ . Below $z\lesssim 1.75$ , there is a high contribution from $\varphi \left(z\right)$ with powers $\alpha =-{2.56}_{-0.54}^{+0.41}$ (broken power law), and $\alpha =-{5.36}_{-2.12}^{+2.60}$ and $\beta =-{2.01}_{-0.50}^{+0.91}$ (triple power law). Meanwhile, for $z\gtrsim 1.75$ , the contribution from $\varphi \left(z\right)$ becomes much less impactful with values $\beta ={0.00}_{-0.31}^{+0.30}$ (broken power law) and $\gamma ={0.06}_{-0.29}^{+0.30}$ (triple power law).

We see the same behavior with the GRB-DD fit. The values for $\beta ,\gamma$ , and $\delta$ that we found coincide with the literature reported values (see Equation (2)) within 1σ (2σ in the case of $\beta$ ). The major difference is in $\alpha$ , with the value found by MCMC being $\alpha ={2.88}_{-1.37}^{+1.28}$ compared to the literature reported value

of $\alpha =0.0157$ (see Equation (2)). Figure 3 shows a plot of the GRB-DD and SFRD for reference. We see that the difference between the two functions is very large at low redshifts but becomes smaller as the redshift increases to $z>4$ .

To summarize, our results show that the Swift data are consistent with the SFRD at high redshifts, and only requires an evolution term that enhances it at low redshifts. This means that if the LGRBs are to follow the SFR, then whatever factors that might affect their detection must have a stronger influence at low redshift than at high redshift. Also, it is notable that high redshift LGRBs detected by Swift are predominantly high luminosity ones. Hence, since our results follow the SFR at high redshifts, they support the findings of [27] .

One important factor to mention that may affect the results is the classification of GRBs. In our analysis, we followed the traditional definition of LGRBs (T₉₀ > 2 s). This definition is widely debated among researchers with many of the newer studies pushing for updated classifications [39] - [45] . Thus, in the future one might investigate this issue further when studying the connection between LGRBs and SFR, by considering a different classification and seeing if it would change the results found in the literature.

4. Conclusion

In this paper, different models for the LGRB redshift distribution were studied. A complete data sample of Swift long GRBs with known redshifts was used without imposing any constraints. Our study shows that the SFRD can describe the LGRB redshift distribution for high redshift but needs an evolution term, to fit the distribution well, at low redshift. Although the data selected may have detection biases, primarily at low redshift, our results support the connection between long GRBs and stellar core-collapse. Our results indicate that the redshift distribution of the data sample fits the SFRD well at high redshifts. However, the issue of excess GRBs at low redshifts persists, because as far as our results go, they do not give us any indication about the nature of the extra contribution term $\varphi \left(z\right)$ . To verify our results, future studies should consider data sets from other satellites or a combined data sample from more than one satellite with a focus on high redshift GRBs.

Appendix A. GRB Dataset

Appendix B. MCMC Contour Plots

NOTES

¹https://swift.gsfc.nasa.gov/archive/grb_table/.

²Note that the best fitting parameters for different cases are expected to be different.

³https://emcee.readthedocs.io/en/stable/.

⁴https://corner.readthedocs.io/en/latest/.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References