A Multilevel Empirical Bayesian Approach to Estimating the Unknown Redshifts of 1366 BATSE Catalog Long-duration Gamma-Ray Bursts

Fundamental to our inference problem is the issue of obtaining a data set of BATSE LGRBs that is minimally biased and representative of the population distribution of LGRBs detectable by the BATSE Large Area Detectors (LADs). The traditional method of GRB classification based on a sharp cutoff line on the observed duration variable ${T}_{90}$ set at $2\mbox{--}3[{\rm{s}}]$ (Kouveliotou et al. 1993) has been shown to be insufficient for an unbiased classification because the duration distributions of LGRBs and SGRBs have significant overlap (e.g., Butler et al. 2010; Shahmoradi 2013a; Shahmoradi & Nemiroff 2015). Instead, we follow the multivariate fuzzy classification approach of Shahmoradi (2013a) and Shahmoradi & Nemiroff (2015) to segregate the two populations of BATSE LGRBs and SGRBs based on their estimated observed spectral peak energies (${E}_{{\rm{p}}}$) from Shahmoradi & Nemiroff (2010) and ${T}_{90}$ from the current BATSE catalog (Goldstein et al. 2013). This leads us to a sample of 1366 LGRBs versus 565 SGRBs in the current BATSE catalog. We refer the interested reader to Shahmoradi (2013a) and Shahmoradi & Nemiroff (2015) for extensive details of the classification procedure.

For our analysis, we also compute, as detailed in Shahmoradi (2013a) and Shahmoradi & Nemiroff (2015), the 1024 [ms] bolometric peak flux (${P}_{\mathrm{bol}}$) and the bolometric fluence (${S}_{\mathrm{bol}}$) of these events from the current BATSE catalog for inclusion in the analysis, in addition to ${E}_{{\rm{p}}}$ and ${T}_{90}$. Taking into account the measurement uncertainties associated with each BATSE event, we can therefore, represent the ith event, ${{\boldsymbol{D}}}_{\mathrm{obs},{\rm{i}}}^{{\rm{g}}}$, in the BATSE catalog by an a priori "known" measurement uncertainty model, ${M}_{{\rm{inois}}}^{{\rm{g}}}$, that together with its "known" parameters, ${M}_{{\mathrm{inois}}^{{\rm{g}}}}^{{\rm{g}}}$, determine the joint four-dimensional probability density function of the observed attributes of the event,

$$\begin{eqnarray}&&\begin{array}{l}\pi \left({{\boldsymbol{D}}}_{\mathrm{obs},{\rm{i}}}^{{\rm{g}}}| {M}_{\mathrm{inois}}^{{\rm{g}}},{M}_{{\mathrm{inois}}^{{\rm{g}}}}^{{\rm{g}}}\right)\\ \quad \propto {M}_{\mathrm{inois}}^{{\rm{g}}}\left({{\boldsymbol{D}}}_{\mathrm{obs},{\rm{i}}}^{{\rm{g}}},{M}_{{\mathrm{inois}}^{{\rm{g}}}}^{{\rm{g}}}\right).\end{array}\end{eqnarray} \tag{ 1 }$$

The modes of these uncertainty models are essentially the values reported in the BATSE catalog and Shahmoradi & Nemiroff (2010),

$$\begin{eqnarray}&&{\widehat{{\boldsymbol{D}}}}_{\mathrm{obs},{\rm{i}}}^{{\rm{g}}}=[{\widehat{P}}_{\mathrm{bol},{\rm{i}}},{\widehat{S}}_{\mathrm{bol},{\rm{i}}},{\widehat{E}}_{{\rm{p}},{\rm{i}}},{\widehat{T}}_{90,{\rm{i}}}],\end{eqnarray} \tag{ 2 }$$

The entire BATSE data set of 1366 LGRB event attributes is then represented by the collection of pairs of such uncertainty models and their parameters,

$$\begin{eqnarray}&&{{ \mathcal D }}_{\mathrm{obs}}^{{\rm{g}}}=\left\{\left({M}_{\mathrm{inois}}^{{\rm{g}}},{M}_{{\mathrm{inois}}^{{\rm{g}}}}^{{\rm{g}}}\right):\,1\leqslant i\leqslant 1366\right\}\end{eqnarray} \tag{ 3 }$$

Note that throughout this manuscript, the appearance of "g" as a superscript solely indicates that the quantity relates to or depends on the four main physical LGRB properties considered in this study, excluding any redshift (z) information.

The BATSE catalog observations appear to have been reported with the assumption that the measurement uncertainty models for all events are multivariate normal distributions, with their mean vectors being the values reported in the catalog and their covariance matrices being diagonal, with the diagonal elements representing the square of the 1–σ errors reported in the catalog.

In sum, our observational data comprises the four main prompt gamma-ray emission properties of LGRBs as described by Equations (1) and (2) and illustrated by Figure 1. ${P}_{\mathrm{bol}}$, as in Equation (2), is, however, a conventionally defined measure of the peak brightness and peripheral to and highly correlated with the more fundamental GRB attribute, ${S}_{\mathrm{bol}}$, its inclusion in our GRB world model is essential as it determines, together with ${E}_{{\rm{p}}}$, the peak photon flux, ${P}_{\mathrm{ph}}$, in 50–300 [keV] range, based upon which BATSE LADs generally triggered on LGRBs.

Zoom In Zoom Out Reset image size

Except for a handful of GRB events, the entire BATSE catalog of GRBs lack any redshift or distance information. Knowledge of individual redshifts is absolutely necessary for accurate cosmographic studies of LGRBs. Essentially, our missing-redshift-data problem for each LGRB reduces to a set of four equations,

$$\begin{eqnarray}\begin{array}{rcl}{\text{}}{L}_{\mathrm{iso}} & = & 4\pi \times {\text{}}{d}_{L}{\left(z\right)}^{2}\times {P}_{\mathrm{bol}},\\ {\text{}}{E}_{\mathrm{iso}} & = & 4\pi \times {\text{}}{d}_{L}{\left(z\right)}^{2}\times {S}_{\mathrm{bol}}/(z+1),\\ {\text{}}{E}_{\mathrm{pz}} & = & {E}_{{\rm{p}}}\times (z+1),\\ {\text{}}{T}_{90{\rm{z}}} & = & {T}_{90}/{\left(z+1\right)}^{\alpha },\end{array}\end{eqnarray} \tag{ 4 }$$

that exactly determine the intrinsic properties of LGRBs: the 1024 [ms] isotropic peak luminosity (${\text{}}{L}_{\mathrm{iso}}$), the total isotropic emission (${\text{}}{E}_{\mathrm{iso}}$), the intrinsic spectral peak energy (${\text{}}{E}_{\mathrm{pz}}$), and the intrinsic duration (${\text{}}{T}_{90{\rm{z}}}$). These four properties are collectively represented by

$$\begin{eqnarray}&&{{\boldsymbol{D}}}_{\mathrm{int},{\rm{i}}}^{{\rm{g}}}=[{L}_{\mathrm{iso},{\rm{i}}},{E}_{\mathrm{iso},{\rm{i}}},{E}_{\mathrm{pz},{\rm{i}}},{T}_{90{\rm{z}},{\rm{i}}}].\end{eqnarray} \tag{ 5 }$$

The term ${\text{}}{d}_{L}(z)$ in Equation (4) represents the luminosity distance,

$$\begin{eqnarray}&&{\text{}}{d}_{L}(z)=\displaystyle \frac{C}{{H}_{0}}(1+z){\int }_{0}^{z}{dz}^{\prime} {\left[{\left(1+z^{\prime} \right)}^{3}{{\rm{\Omega }}}_{M}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}\right]}^{-1/2},\end{eqnarray} \tag{ 6 }$$

and the exponent $\alpha =0.66$ in the mapping of ${T}_{90}$ to ${\text{}}{T}_{90{\rm{z}}}$ takes into account the cosmological time dilation as well as an energy-band correction (i.e., K-correction) of the form ${\left(1+z\right)}^{-0.34}$ to the observed durations (e.g., Gehrels et al. 2006). Throughout this work, we assume a flat ΛCDM cosmology, with parameters set to h = 0.70, ${{\rm{\Omega }}}_{}M=0.27$, and ${{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.73$ (Jarosik et al. 2011). Also, the parameters C and ${H}_{0}=100h$ [km s^–1 MPc^–1] stand for the speed of light and the Hubble constant, respectively. Equations (4) can be linearized by taking the logarithms of both sides, leading to

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}({\text{}}{L}_{\mathrm{iso}}) & = & \mathrm{log}({P}_{\mathrm{bol}})+2\mathrm{log}({\text{}}{d}_{L})+\mathrm{log}(4\pi ),\\ \mathrm{log}({\text{}}{E}_{\mathrm{iso}}) & = & \mathrm{log}({S}_{\mathrm{bol}})+2\mathrm{log}({\text{}}{d}_{L})+\mathrm{log}(4\pi )-\mathrm{log}(z+1),\\ \mathrm{log}({\text{}}{E}_{\mathrm{pz}}) & = & \mathrm{log}({E}_{{\rm{p}}})+\mathrm{log}(z+1),\\ \mathrm{log}({\text{}}{T}_{90{\rm{z}}}) & = & \mathrm{log}({T}_{90})-\alpha \mathrm{log}(z+1).\end{array}\end{eqnarray} \tag{ 7 }$$

More concisely, we can write the four Equations (7) for the ith GRB event as a single equation,

$$\begin{eqnarray}&&\mathrm{log}({{\boldsymbol{D}}}_{\mathrm{int},{\rm{i}}}^{{\rm{g}}})=\mathrm{log}({{\boldsymbol{D}}}_{\mathrm{obs},{\rm{i}}}^{{\rm{g}}})+{ \mathcal F }({z}_{i}),\end{eqnarray} \tag{ 8 }$$

where ${ \mathcal F }({z}_{i})$ describes the mapping from the observer to the rest frame of the GRB given its known redshift, z_i. This simple equation deserves some deliberations. Despite its simplicity, this single equation is algebraically unsolvable for almost any BATSE catalog GRB as it contains two unknowns (${{\boldsymbol{D}}}_{\mathrm{int},{\rm{i}}}^{{\rm{g}}}$, z_i).

Although Equation (8) is algebraically degenerate, the two unknowns of the equation can still be constrained probabilistically. To understand how this is possible, consider the following simple toy problem.

Without loss of generality, suppose the observed properties of individual BATSE LGRBs are exactly known, with no measurement error, as illustrated by the individual blue-colored vertical lines in the bottom plot of Figure 2. The corresponding redshifts of these events, (represented by the black lines in the middle plot of Figure 2) are, however, unknown, and we wish to estimate them. Although we have no knowledge of the joint population distribution of the intrinsic properties of LGRBs, illustrated by the red distribution on the top plot of Figure 2, there are strong arguments in favor of these properties being potentially well described by a four-dimensional multivariate log-normal distribution, ${ \mathcal N }({\boldsymbol{\mu }},{\boldsymbol{\Sigma }})$, in the space of ${\text{}}{L}_{\mathrm{iso}}$, ${\text{}}{E}_{\mathrm{iso}}$, ${\text{}}{E}_{\mathrm{pz}}$, and ${\text{}}{T}_{90{\rm{z}}}$ (Shahmoradi 2013a; Shahmoradi & Nemiroff 2015). Here, ${\boldsymbol{\mu }}$ and ${\boldsymbol{\Sigma }}$ represent the mean and the covariance matrix of the multivariate log-normal distribution, respectively.

Zoom In Zoom Out Reset image size

We now reach the crucial step in the inference process: despite the complete lack of information about the redshifts of BATSE LGRBs, we can use the existing prior knowledge about the overall cosmic redshift distribution of LGRBs to integrate over all possible redshifts for each observed LGRB in the BATSE catalog to infer a range of plausible values for the intrinsic properties of the corresponding LGRB. These individually computed probability density functions (PDFs) of the intrinsic properties can be then used to infer the unknown parameters $({\boldsymbol{\mu }},{\boldsymbol{\Sigma }})$ of the joint population distribution of the intrinsic properties of LGRBs (i.e., the multivariate log-normal distribution).

Once $({\boldsymbol{\mu }},{\boldsymbol{\Sigma }})$ are constrained, we can use the inferred population distribution of the intrinsic LGRB properties together with the observed properties to estimate the redshifts of individual BATSE LGRBs, independently of each other. The estimated redshifts can again be used to further constrain ${ \mathcal N }({\boldsymbol{\mu }},{\boldsymbol{\Sigma }})$, which will then result in even tighter estimates for the individual redshifts of BATSE LGRBs. This recursive progress can practically continue until convergence to a set of fixed individual redshift estimates occurs.

At first glance, this simple semi-Bayesian mathematical approach may sound like magic and perhaps, too good to be true. Sometimes it is. However, as we later explain in Section 4, it can also lead to reasonably accurate results if some conditions regarding the problem and the observational data set are satisfied.

More formally, we model the process of LGRB observation as a nonhomogeneous Poisson process whose mean rate parameter is the "censored" cosmic LGRB rate, ${{ \mathcal R }}_{\mathrm{cen}}$. Representing each LGRB by

$$\begin{eqnarray}&&{{\boldsymbol{D}}}_{\mathrm{int},{\rm{i}}}=\{{{\boldsymbol{D}}}_{\mathrm{int},{\rm{i}}}^{{\rm{g}}},{z}_{{\rm{i}}}\},\,1\leqslant i\leqslant 1366,\end{eqnarray} \tag{ 9 }$$

where ${{\boldsymbol{D}}}_{\mathrm{int},{\rm{i}}}^{{\rm{g}}}$ is defined by Equation (5), we compute the probability of occurrence of each BATSE LGRB event in the five-dimensional attributes space, ${\rm{\Omega }}({{\boldsymbol{D}}}_{\mathrm{int}})$, of z, ${\text{}}{L}_{\mathrm{iso}}$, ${\text{}}{E}_{\mathrm{iso}}$, ${\text{}}{E}_{\mathrm{pz}}$, ${\text{}}{T}_{90{\rm{z}}}$, as a function of the parameters, ${{\boldsymbol{\theta }}}_{\mathrm{cen}}$, of the observed LGRB rate model, ${{ \mathcal R }}_{\mathrm{cen}}$,

$$\begin{eqnarray}&&\pi \left({{\boldsymbol{D}}}_{\mathrm{int},{\rm{i}}}| {{ \mathcal R }}_{\mathrm{cen}},{{\boldsymbol{\theta }}}_{\mathrm{cen}}\right)\propto {{ \mathcal R }}_{\mathrm{cen}}\left({{\boldsymbol{D}}}_{\mathrm{int},{\rm{i}}},{{\boldsymbol{\theta }}}_{\mathrm{cen}}\right),\end{eqnarray} \tag{ 10 }$$

where ${{ \mathcal R }}_{\mathrm{cen}}$ represents the BATSE-censored rate of LGRB occurrence in the universe (due to the BATSE detection efficiency limitations as detailed in Section 2.6). This equation can be further expanded in terms of BATSE detection efficiency function, ${\eta }_{\mathrm{eff}}$, and the true cosmic LGRB rate, ${{ \mathcal R }}_{\mathrm{tru}}$, as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dN}}_{\mathrm{obs}}}{d{{\boldsymbol{D}}}_{\mathrm{int}}} & = & {{ \mathcal R }}_{\mathrm{cen}}\left({{\boldsymbol{D}}}_{\mathrm{int}},{{\boldsymbol{\theta }}}_{\mathrm{cen}}\right),\\ & = & {\eta }_{\mathrm{eff}}\left({{\boldsymbol{D}}}_{\mathrm{int}},{{\boldsymbol{\theta }}}_{\mathrm{eff}}\right)\times {{ \mathcal R }}_{\mathrm{tru}}\left({{\boldsymbol{D}}}_{\mathrm{int}},{{\boldsymbol{\theta }}}_{\mathrm{tru}}\right)\end{array}\end{eqnarray} \tag{ 11 }$$

for a given set of input intrinsic LGRB attributes, ${{\boldsymbol{D}}}_{\mathrm{int}}$, with ${{\boldsymbol{\theta }}}_{\mathrm{cen}}=\{{{\boldsymbol{\theta }}}_{\mathrm{eff}},{{\boldsymbol{\theta }}}_{\mathrm{tru}}\}$ as the set of the parameters of our models for the BATSE detection efficiency and the intrinsic cosmic LGRB rate, respectively. Assuming no systematic evolution of LGRB characteristics with redshift, which is a plausible assumption supported by independent studies (e.g., Butler et al. 2010; hereafter B10), the intrinsic LGRB rate itself can be written as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{{dN}}_{\mathrm{int}}}{d{{\boldsymbol{D}}}_{\mathrm{int}}} & = & {{ \mathcal R }}_{\mathrm{tru}}\left({{\boldsymbol{D}}}_{\mathrm{int}},{{\boldsymbol{\theta }}}_{\mathrm{tru}}\right)\\ & = & {{ \mathcal R }}_{\mathrm{tru}}^{{\rm{g}}}\left({{\boldsymbol{D}}}_{\mathrm{int}}^{{\rm{g}}},{{\boldsymbol{\theta }}}_{\mathrm{tru}}^{{\rm{g}}}\right)\times \displaystyle \frac{\dot{\zeta }(z,{{\boldsymbol{\theta }}}_{z}){dV}/{dz}}{(1+z)},\end{array}\end{eqnarray} \tag{ 12 }$$

with ${{\boldsymbol{\theta }}}_{\mathrm{tru}}=\{{{\boldsymbol{\theta }}}_{\mathrm{tru}}^{{\rm{g}}},{{\boldsymbol{\theta }}}_{z}\}$, where ${{ \mathcal R }}_{\mathrm{tru}}^{{\rm{g}}}$ is a statistical model, with ${{\boldsymbol{\theta }}}_{\mathrm{tru}}^{{\rm{g}}}$ denoting its parameters that describe the population distribution of LGRBs in the four-dimensional attributes space of ${{\boldsymbol{D}}}_{\mathrm{int}}^{{\rm{g}}}=[{\text{}}{L}_{\mathrm{iso}},{\text{}}{E}_{\mathrm{iso}},{\text{}}{E}_{\mathrm{pz}},{\text{}}{T}_{90{\rm{z}}}]$, and the term $\dot{\zeta }(z,{{\boldsymbol{\theta }}}_{z})$ represents the comoving rate density model of LGRBs with the set of parameters ${{\boldsymbol{\theta }}}_{z}$, while the factor $(1+z)$ in the denominator accounts for the cosmological time dilation. The comoving volume element per unit redshift, ${dV}/{dz}$, is given by (e.g., Winberg 1972; Peebles 1993)

$$\begin{eqnarray}&&\displaystyle \frac{{dV}}{{dz}}=\displaystyle \frac{C}{{H}_{0}}\displaystyle \frac{4\pi {\text{}}{d}_{L}^{2}(z)}{{\left(1+z\right)}^{2}{\left[{{\rm{\Omega }}}_{{\rm{M}}}{\left(1+z\right)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}\right]}^{1/2}},\end{eqnarray} \tag{ 13 }$$

with ${\text{}}{d}_{L}$ standing for the luminosity distance as given in Equation (6). If the three rate models, $(\dot{\zeta },{\eta }_{\mathrm{eff}},{{ \mathcal R }}_{\mathrm{tru}}^{{\rm{g}}})$, and their parameters were known a priori, one could readily compute the PDFs of the set of unknown redshifts of all BATSE LGRBs,

$$\begin{eqnarray}&&{ \mathcal Z }=\left\{{z}_{{\rm{i}}}:\,1\leqslant i\leqslant 1366\right\},\end{eqnarray} \tag{ 14 }$$

as

$$\begin{eqnarray}&&\begin{array}{l}\pi \left({ \mathcal Z }| {{ \mathcal D }}_{\mathrm{obs}}^{{\rm{g}}},{{ \mathcal R }}_{\mathrm{cen}},{{\boldsymbol{\theta }}}_{\mathrm{cen}}\right)\\ \quad \propto {\int }_{{\rm{\Omega }}({{ \mathcal D }}_{\mathrm{int}}^{{\rm{g}}})}{{ \mathcal R }}_{\mathrm{cen}}\left({ \mathcal Z },{{{ \mathcal D }}_{\mathrm{int}}^{{\rm{g}}}}^{* },{{\boldsymbol{\theta }}}_{\mathrm{cen}}\right)\,d{{{ \mathcal D }}_{\mathrm{int}}^{{\rm{g}}}}^{* },\end{array}\end{eqnarray} \tag{ 15 }$$

where the integration is performed over all possible realizations, ${{{ \mathcal D }}_{\mathrm{int}}^{{\rm{g}}}}^{* }$, of the BATSE LGRB data set given the measurement uncertainty models in Equation (3). For a range of possible parameter values, the redshift probabilities can be computed by marginalizing over the entire parameter space, ${\rm{\Omega }}({{\boldsymbol{\theta }}}_{\mathrm{cen}})$, of the model,

$$\begin{eqnarray}&&\begin{array}{l}\pi \left({ \mathcal Z }| {{ \mathcal D }}_{\mathrm{obs}}^{{\rm{g}}},{{ \mathcal R }}_{\mathrm{cen}}\right)\\ ={\int }_{{\rm{\Omega }}({{\boldsymbol{\theta }}}_{\mathrm{cen}})}\pi \left({ \mathcal Z }| {{ \mathcal D }}_{\mathrm{obs}}^{{\rm{g}}},{{ \mathcal R }}_{\mathrm{cen}},{{\boldsymbol{\theta }}}_{\mathrm{cen}}\right)\\ \quad \times \,\pi \left({{\boldsymbol{\theta }}}_{\mathrm{cen}}| {{ \mathcal D }}_{\mathrm{obs}}^{{\rm{g}}},{{ \mathcal R }}_{\mathrm{cen}}\right)\,d{{\boldsymbol{\theta }}}_{\mathrm{cen}}.\end{array}\end{eqnarray} \tag{ 16 }$$

The problem, however, is that neither the rate models nor their parameters are known a priori. Even more problematic is the circular dependency of the posterior PDFs of ${ \mathcal Z }$ and ${{\boldsymbol{\theta }}}_{\mathrm{cen}}$ on each other,

$$\begin{eqnarray}&&\begin{array}{l}\pi \left({{\boldsymbol{\theta }}}_{\mathrm{cen}}\,| \,{{ \mathcal D }}_{\mathrm{obs}}^{{\rm{g}}},{{ \mathcal R }}_{\mathrm{cen}}\right)\\ ={\int }_{{\rm{\Omega }}({ \mathcal Z })}\pi \left({{\boldsymbol{\theta }}}_{\mathrm{cen}}| { \mathcal Z },{{ \mathcal D }}_{\mathrm{obs}}^{{\rm{g}}},{{ \mathcal R }}_{\mathrm{cen}}\right)\\ \quad \times \,\pi \left({ \mathcal Z }\,| \,{{ \mathcal D }}_{\mathrm{obs}}^{{\rm{g}}},{{ \mathcal R }}_{\mathrm{cen}}\right)\,d{ \mathcal Z }.\end{array}\end{eqnarray} \tag{ 17 }$$

Therefore, we adopt the following methodology, which is reminiscent of the Empirical Bayes (Robbins 1985) and Expectation-Maximization algorithms (Dempster et al. 1977), to estimate the redshifts of BATSE LGRBs. First, we propose models for $(\dot{\zeta },{\eta }_{\mathrm{eff}},{{ \mathcal R }}_{\mathrm{tru}}^{{\rm{g}}})$, whose parameters have yet to be constrained by observational data. Given the three rate models, we can then proceed to constrain the free parameters of the observed cosmic LGRB rate, ${{ \mathcal R }}_{\mathrm{cen}}$, based on BATSE LGRB data.

The most appropriate fitting approach should take into account the observational uncertainties and any prior knowledge from independent sources. This can be achieved via the multilevel Bayesian methodology (e.g., Shahmoradi 2017) by constructing the likelihood function and the posterior PDF of the parameters of the model, while taking into account the uncertainties in observational data (e.g., Equation (61) in Shahmoradi 2017),

$$\begin{eqnarray}&&\begin{array}{l}\pi \left({{\boldsymbol{\theta }}}_{\mathrm{cen}}| {{ \mathcal D }}_{\mathrm{obs}}^{{\rm{g}}},{{ \mathcal R }}_{\mathrm{cen}}\right)=\displaystyle \frac{\pi \left({{\boldsymbol{\theta }}}_{\mathrm{cen}}| {{ \mathcal R }}_{\mathrm{cen}}\right)}{\pi \left({{ \mathcal D }}_{\mathrm{obs}}^{{\rm{g}}}| {{ \mathcal R }}_{\mathrm{cen}}\right)}\\ \,{\int }_{{\rm{\Omega }}({ \mathcal Z })}{\int }_{{\rm{\Omega }}({{ \mathcal D }}_{\mathrm{int}}^{{\rm{g}}})}\pi \left({{{ \mathcal D }}_{\mathrm{int}}^{{\rm{g}}}}^{* }\,| {{ \mathcal D }}_{\mathrm{obs}}^{{\rm{g}}},{ \mathcal Z },{{ \mathcal R }}_{\mathrm{cen}},{{\boldsymbol{\theta }}}_{\mathrm{cen}}\right)\pi \left({ \mathcal Z }| {{ \mathcal R }}_{\mathrm{cen}},{{\boldsymbol{\theta }}}_{\mathrm{cen}}\right)\\ \quad \times \,d{{{ \mathcal D }}_{\mathrm{int}}^{{\rm{g}}}}^{* }\,d{ \mathcal Z },\end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\begin{array}{l}\mathop{\propto }\limits^{{\rm{i}}.{\rm{i}}.{\rm{d}}.}\exp \left(-{\int }_{{\rm{\Omega }}({{\boldsymbol{D}}}_{\mathrm{int}})}{{ \mathcal R }}_{\mathrm{cen}}\left({{\boldsymbol{D}}}_{\mathrm{int}}^{* },{{\boldsymbol{\theta }}}_{\mathrm{cen}}\right)d{{\boldsymbol{D}}}_{\mathrm{int}}^{* }\right)\\ \quad \times \,\displaystyle \prod _{i=1}^{1366}{\eta }_{\mathrm{eff}}\left({\widehat{{\boldsymbol{D}}}}_{\mathrm{obs},{\rm{i}}}^{{\rm{g}}},{{\boldsymbol{\theta }}}_{\mathrm{eff}}\right){\int }_{{\rm{\Omega }}({{\boldsymbol{D}}}_{\mathrm{int}})}{{ \mathcal R }}_{\mathrm{tru}}\left({{\boldsymbol{D}}}_{\mathrm{int}}^{* },{{\boldsymbol{\theta }}}_{\mathrm{tru}}\right)\\ \quad \times \,\pi \left({{\boldsymbol{D}}}_{\mathrm{int}}^{* }| \,{{\boldsymbol{D}}}_{\mathrm{obs},{\rm{i}}}^{{\rm{g}}}\right)d{{\boldsymbol{D}}}_{\mathrm{int}}^{* },\end{array}\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\begin{array}{l}\mathop{\equiv }\limits^{\mathrm{no\; noise}}\exp \left(-{\int }_{{\rm{\Omega }}({{\boldsymbol{D}}}_{\mathrm{int}})}{{ \mathcal R }}_{\mathrm{cen}}\left({{\boldsymbol{D}}}_{\mathrm{int}}^{* },{{\boldsymbol{\theta }}}_{\mathrm{cen}}\right)d{{\boldsymbol{D}}}_{\mathrm{int}}^{* }\right)\\ \times \,\displaystyle \prod _{i=1}^{1366}{\eta }_{\mathrm{eff}}\left({\widehat{{\boldsymbol{D}}}}_{\mathrm{obs},{\rm{i}}}^{{\rm{g}}},{{\boldsymbol{\theta }}}_{\mathrm{eff}}\right){\int }_{{z}^{* }=0}^{{z}^{* }=+\infty }{{ \mathcal R }}_{\mathrm{tru}}\left({\widehat{{\boldsymbol{D}}}}_{\mathrm{obs},{\rm{i}}}^{{\rm{g}}},{z}^{* },{{\boldsymbol{\theta }}}_{\mathrm{tru}}\right){{dz}}^{* },\end{array}\end{eqnarray} \tag{ 20 }$$

where Equation (19) holds under the assumption of an independent and identical distribution (i.e., the i.i.d. property) of BATSE LGRBs, and the second integration within it is performed over all possible realizations of the truth for the ith observed BATSE LGRB, ${{\boldsymbol{D}}}_{\mathrm{obs},{\rm{i}}}^{{\rm{g}}}$. Equation (19) can be further considerably simplified to Equation (20) by assuming no measurement uncertainty in the observational data, except redshift (z), which is completely unknown for BATSE LGRBs.

Once the posterior PDF of the model parameters is obtained, it can be plugged into Equation (16) to constrain the redshift PDF of individual BATSE LGRBs at the second level of modeling.

Our main assumption in this work is that the intrinsic comoving rate density of LGRBs closely traces the comoving SFR density (e.g., Madau & Dickinson 2014; Madau & Fragos 2017; Fermi-LAT Collaboration et al. 2018; hereafter M14, M17, F18, respectively) or a metallicity-corrected SFR density as prescribed by B10. Consequently, regardless of the individually unknown redshifts of BATSE LGRBs, the overall redshift distribution of all BATSE LGRBs together is enforced in our modeling to follow the cosmic SFR convolved with the BATSE detection efficiency model, ${\eta }_{\mathrm{eff}}$, which is detailed in Section 2.6. This assumption is essential for the success of our modeling approach, as any attempts to constrain the comoving rate density, $\dot{\zeta }(z)$, of LGRBs solely based on BATSE data leads to highly degenerate parameter space, ${\rm{\Omega }}({{\boldsymbol{\theta }}}_{\mathrm{cen}})$, and parameter estimates for our model, ${{ \mathcal R }}_{\mathrm{cen}}$.

As for the choice of the LGRB rate density model, $\dot{\zeta }$, we have considered and simulated six different LGRB rate density scenarios, three of which have the generic continuous piecewise form,

$$\begin{eqnarray}\dot{\zeta }(z)\propto \left\{\begin{array}{ll}{\left(1+z\right)}^{{\gamma }_{0}} & z\lt {z}_{0}\\ {\left(1+z\right)}^{{\gamma }_{1}} & {z}_{0}\lt z\lt {z}_{1}\\ {\left(1+z\right)}^{{\gamma }_{2}} & z\gt {z}_{1},\end{array}\right.\end{eqnarray} \tag{ 21 }$$

with parameters

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{\theta }}}_{z} & = & ({z}_{0},{z}_{1},{\gamma }_{0},{\gamma }_{1},{\gamma }_{2})\\ & = & \left\{\begin{array}{ll}(0.97,4.5,3.4,-0.3,-7.8) & ({\rm{H}}06)\\ (0.993,3.8,3.3,0.055,-4.46) & ({\rm{L}}08)\\ (0.97,4.00,3.14,1.36,-2.92) & ({\rm{B}}10)\\ (0.0,4.5,0.0,0.0,-7.8) & ({\rm{P}}15)\end{array}\right.\end{array}\end{eqnarray} \tag{ 22 }$$

corresponding to the SFR density of Hopkins & Beacom (2006, hereafter H06), Li (2008, hereafter L08), and Petrosian et al. (2015, hereafter P15), and a bias-corrected redshift distribution of LGRBs derived from Swift data by B10. The other three redshift scenarios follow the generic functional of Equation (15) in Madau & Dickinson (2014),

$$\begin{eqnarray}&&\dot{\zeta }(z)\propto \displaystyle \frac{{\left(1+z\right)}^{{\alpha }_{1}}}{{\left(1+z\right)}^{{\alpha }_{2}}+{\left(1+{z}_{0}\right)}^{{\alpha }_{2}}},\end{eqnarray} \tag{ 23 }$$

with parameters

$$\begin{eqnarray}{{\boldsymbol{\theta }}}_{z}=({z}_{0},{\alpha }_{1},{\alpha }_{2})=\left\{\begin{array}{ll}(1.90,2.70,5.60) & ({\rm{M}}14)\\ (2.20,2.60,6.20) & ({\rm{M}}17)\\ (1.63,2.99,6.19) & ({\rm{F}}18)\end{array}\right..\end{eqnarray} \tag{ 24 }$$

For the sake of brevity, we will only present the results for four out of these seven SFR density scenarios: H06, B10, P15, M17.

We model the joint four-dimensional distribution of ${{\boldsymbol{D}}}_{\mathrm{int}}^{{\rm{g}}}$ with a multivariate log-normal distribution, ${{ \mathcal R }}_{\mathrm{tru}}^{{\rm{g}}}\equiv { \mathcal L }{ \mathcal N }$, whose parameters (i.e., the mean vector and the covariance matrix), ${{\boldsymbol{\theta }}}_{\mathrm{tru}}^{{\rm{g}}}=\{{\boldsymbol{\mu }},{\boldsymbol{\Sigma }}\}$, will have to be constrained by data. The justification for the choice of a multivariate log-normal as the underlying intrinsic population distribution of LGRBs is multifold. First, the observed joint distribution of BATSE LGRB properties highly resembles a log-normal shape that is censored close to the detection threshold of BATSE. Second, unlike the power-law distribution which has traditionally been the default choice of model for the luminosity function of LGRBs, log-normal models provide natural upper and lower bounds on the total energy budget and luminosity of LGRBs, eliminating the need for setting artificially sharp bounds on the distributions to properly normalize them. Third, a log-normal along with a Gaussian distribution is among the most naturally occurring statistical distributions in nature, whose generalizations to multiple dimensions are also well studied and understood. This is a highly desired property especially for our work, given the overall mathematical and computational complexity of the model proposed and developed here.

Compared to the Fermi Gamma-ray Burst Monitor (Meegan et al. 2009) and Neil Gehrels Swift Observatory (Gehrels et al. 2004; Lien et al. 2016), BATSE had a relatively simple triggering algorithm. The BATSE detection efficiency and algorithm has been already extensively studied by the BATSE team as well as independent authors (e.g., Pendleton et al. 1998, 1995; Hakkila et al. 2003; see Shahmoradi & Nemiroff 2010, 2011a; Shahmoradi 2013a; Shahmoradi & Nemiroff 2015 for further discussion and references). However, simple implementation and usage of the known BATSE trigger threshold for modeling the BATSE catalog's sample incompleteness can lead to systematic biases in the inferred quantities of interest. BATSE triggered on 2702 GRBs, out of which only 2145, or approximately 79%, have been consistently analyzed and reported in the current BATSE catalog, with the remaining 21% either having a low accumulation of count rates or missing full spectral/temporal coverage (Goldstein et al. 2013). Thus, the extent of sample incompleteness in the BATSE catalog is likely not fully and accurately represented by the BATSE triggering algorithm alone.

BATSE LADs generally triggered on a GRB if the number of photons per 1024 [ms] arriving at the detectors in the 50–300 [keV] energy window, ${P}_{\mathrm{ph}}$, reached a certain threshold in units of the background photon count fluctuations, σ. This threshold was typically set to $5.5\sigma$ during much of BATSE's operational lifetime. However, the naturally occurring fluctuations in the average background photon counts effectively lead to a monotonically increasing BATSE detection efficiency as a function of ${P}_{\mathrm{ph}}$ instead of a sharp cutoff on the observed ${P}_{\mathrm{ph}}$ distribution of LGRBs. Therefore, we model the effects of BATSE detection efficiency and sample incompleteness more accurately by an error function,

$$\begin{eqnarray}&&\begin{array}{l}\pi \left(\mathrm{detection}| {\mu }_{\mathrm{th}},{\sigma }_{\mathrm{th}},{P}_{\mathrm{ph}}\right)\\ \quad =\,\displaystyle \frac{1}{2}+\displaystyle \frac{1}{2}\times \mathrm{erf}\left(\displaystyle \frac{{\mathrm{log}}_{10}{P}_{\mathrm{ph}}-{\mu }_{\mathrm{th}}}{{\sigma }_{\mathrm{th}}\sqrt{2}}\right),\end{array}\end{eqnarray} \tag{ 25 }$$

that, for a given set of the error function's parameters (${{\boldsymbol{\theta }}}_{\mathrm{eff}}=\{{\mu }_{\mathrm{th}},{\sigma }_{\mathrm{th}}\}$), yields the probability of the detection of an LGRB with ${P}_{\mathrm{ph}}$ photon counts per second. Due to the unknown effects of sample incompleteness on BATSE data, we leave these two threshold parameters free solely to be constrained by the observational data. Figure 3 compares the resulting best-fit detection efficiency functions for BATSE under different SFR scenarios with the BATSE nominal detection efficiency for LGRBs⁴ (see also Pendleton et al. 1995, 1998; Paciesas et al. 1999; Hakkila et al. 2003; Shahmoradi 2013a; Shahmoradi & Nemiroff 2015).

Zoom In Zoom Out Reset image size

The current BATSE catalog already provides estimates of ${P}_{\mathrm{ph}}$ for all 1366 LGRBs in our analysis. The connection between ${P}_{\mathrm{ph}}$ and the bolometric 1024 [ms] peak flux, ${P}_{\mathrm{bol}}$, which is used in our modeling, is provided by fitting all LGRB spectra with a smoothly broken power law,

$$\begin{eqnarray}{\rm{\Phi }}(E)\propto \left\{\begin{array}{ll}{E}^{\alpha }{{\rm{e}}}^{\left(-\tfrac{(1+z)(2+\alpha )E}{{\text{}}{E}_{\mathrm{pz}}}\right)} & \mathrm{if}\,E\leqslant \left(\displaystyle \frac{{\text{}}{E}_{\mathrm{pz}}}{1+z}\right)\left(\tfrac{\alpha -\beta }{2+\alpha }\right)\\ {E}^{\beta } & \mathrm{if}\,\mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 26 }$$

known as the Band model (Band et al. 1993) to infer their spectral normalization constants while fixing the low- and high-energy photon indices of the Band model to the corresponding population averages $(\alpha ,\beta )=(-1.1,-2.3)$ and fixing the observed spectral peak energies, ${E}_{{\rm{p}}}$, of individual bursts to the corresponding best-fit values from Shahmoradi & Nemiroff (2010). Such an approximation is reasonable given the typically large uncertainties that exist in the spectral and temporal properties of GRBs (e.g., Butler et al. 2010; Shahmoradi 2013a) and the relatively large variance of the population distribution of ${P}_{\mathrm{bol}}$ in our sample.

A remarkable modeling assumption in our work is the redshift independence of the proposed model for the four main gamma-ray properties of cosmic LGRBs (the term ${{ \mathcal R }}_{\mathrm{tru}}^{{\rm{g}}}$ in Equation (12)). This modeling assumption is similar to the assumptions of Shahmoradi (2013a) and Shahmoradi & Nemiroff (2015) and follows the findings of B10, who, based on a multivariate analysis of a large sample of Swift LGRBs with and without redshifts, reject the hypothesis of an evolving redshift-dependent luminosity function for LGRBs. This is in contrast to the findings of Petrosian et al. (2015), Yu et al. (2015), Pescalli et al. (2016), Tsvetkova et al. (2017), and Lloyd-Ronning et al. (2019), who model the cosmic rates of LGRBs via a redshift-dependent luminosity function.

Regardless of the validity of the aforementioned assumption in our modeling, we note that it is practically impossible to infer any redshift dependence of the energetics and the luminosity function of LGRBs solely via BATSE LGRB data. This limitation primarily results from the complete lack of knowledge of the redshifts of BATSE LGRBs. This missing knowledge leads to a highly degenerate parameter space for models that consider the possibility of a redshift evolution of the prompt emission properties of LGRBs.

The multivariate log-normal luminosity function assumption of this work also needs to be further studied and validated, and its effects on the inferred redshifts understood better in a future work. As illustrated in Figure 2, the accuracy of the inferred redshifts of LGRBs depends strongly on the width of the redshift distribution relative to the intrinsic distributions of LGRB properties, most importantly, the ${\text{}}{L}_{\mathrm{iso}}$ and ${\text{}}{E}_{\mathrm{iso}}$ distributions. If the intrinsic LGRB property distributions are too broad compared with the redshift distribution, then one would expect less accurate redshift estimates. However, we anticipate that the particular choice of the luminosity function (for example, log-normal versus power law or broken power law) will have minimal effects on the inferred redshifts because the detection threshold of gamma-ray telescopes cuts the intrinsic LGRB property distributions before any noticeable statistically significant discrepancies between the different choices of energetics distributions could be identified solely from the observational data. This is an important issue that needs further investigation and validation in the future.

The investigation of the evidence for the luminosity evolution of LGRBs goes beyond the scope of this manuscript. We instead defer a detailed investigation of this problem to a future paper. Nevertheless, we emphasize that the results presented in this manuscript are completely based on the validity of the assumption of the redshift independence of the luminosity function of LGRBs (e.g., Butler et al. 2010) that has been recently challenged by independent studies (e.g., Petrosian et al. 2015; Yu et al. 2015; Pescalli et al. 2016; Tsvetkova et al. 2017; Lloyd-Ronning et al. 2019).