© ESO 2022

1. Introduction

Gamma-ray bursts (GRBs) are among the most energetic objects in the Universe. After their discovery (Klebesadel et al. 1973), it was shown that there are two types of GRBs (Mazets et al. 1981): short and hard, and long and soft. This is well supported by both observational data and statistical analysis (for a survey and the relevant references see, e.g., Mészáros 2006).

An intrinsic luminosity function (ILF) is among the most sought-after quantities for any class of astrophysical objects. The intrinsic luminosity function (ILF) of GRBs has been the subject of many papers in the past, but this subject has been plagued by ignorance of what its physical significance would be even if it were known. In particular, as GRBs are known to be beamed along a jet axis, the intrinsic luminosity function also sheds light on the angular emission profile of the GRBs, assuming the observers are distributed isotropically.

There exists inevitable scatter in the intrinsic luminosities of GRBs, and it is not clear whether or not this scatter is due to relativistic beaming, as offset viewers may see diminished luminosity. Scatter in the intrinsic luminosity can cause selection bias. When the scatter in luminosity is strong, it is important to quantify the effect in order to distinguish it from true dynamical evolution. Therefore, it is important to study luminosity selection effects in GRBs in order to identify the bias.

Several attempts have been made in the past to find the ILF by fitting the peak flux distribution and the redshift distribution (Tsvetkova et al. 2017; Petrosian et al. 2015; Pescalli et al. 2016; Guetta & Piran 2006). A common result from all these studies is that the ILF that best fits the data is a broken power law with three free parameters (the low and high power-law index and the luminosity break). Guetta et al. (2005) also tried to get information on the angular distribution of the GRBs by assuming that, in the case of a uniform jet, the quantity Lθ² may be considered constant. This assumption was implied by the observation on the sample considered by Bloom et al. (2003).

The question, related to finding the ILF, has even greater significance when considering the ILF of ‘low-luminosity’ GRBs (i.e. L ≤ 10⁴⁹ erg s⁻¹). Are they low luminosity because they are viewed from far away from the direction of motion of – and hence the direction of beaming by – the material with which the observed photons last interacted, or because material beaming mostly in the direction of the observer is less luminous? Guetta & Della Valle (2007) find that there is a difference in the frequency of occurrence of high-luminosity GRBs and low-luminosity GRBs. This could be interpreted as the existence of two physically distinct classes of GRBs (Cobb et al. 2006; Soderberg et al. 2006) and whereby low-luminosity GRBs are (intrinsically) more frequent events than high-luminosity GRBs. Alternatively, there could simply be a single population whose rate is luminosity dependent. This single population originates in both the isotropic and subenergetic component, detectable only in nearby GRBs, and the highly collimated component, observable by sampling huge volumes of space and thus mainly detectable in high-z GRBs.

In this paper, we propose a physical hypothesis for the ILF. We test this physical hypothesis against observed data using statistical tools based on median luminosity where the best-fit model is the one that divides the data most symmetrically.

The diagnostic we suggest is based on comparison of the median luminosity L_median of a set of GRBs that have a known intrinsic luminosity function with the observed luminosity of a wide range of GRBs. The advantage of this method is that it requires only a small data set.

Though we perform the statistical analysis for a large GRB sample covering both low and high redshifts, the diagnostic proposed in this paper is well suited even for small sets of GRBs, such as low-redshift GRBs, z ≤ 0.15, of which only about 10 have been discovered, as long as the set is not chosen with luminosity bias.

Several authors have shown that the luminosity function evolves significantly with redshift (Petrosian et al. 2015; Pescalli et al. 2016; Lloyd-Ronning et al. 2019). However, the method proposed in this paper is independent of the redshift distribution. Moreover, we show that it is also independent of the maximum luminosity of GRBs. The paper is organized as follows: In Sect. 2, we discuss our choice of mathematical form of luminosity function, and its physical motivation. In Sect. 3, we outline how to calculate the median bolometric luminosity L_median of detected sources as a function of redshift given some intrinsic luminosity for all of them. In Sect. 3, we compare L_median calculated from a set of possibilities, and show that the data are consistent with our physical hypotheses. Finally, in Sect. 4 we report our conclusion.

2. The physical model and implications for the form of the GRB ILF

We consider a model of a GRB jet with constant intrinsic luminosity, beaming angle θ₀, and bulk Lorentz factor Γ. We are aware that these assumptions are questionable. There have been many papers (Shuang-Xi et al. 2017; Ackermann et al. 2011; Perna & Loeb 1998; Zou et al. 2017; Ghirlanda et al. 2013; Tang et al. 2015) dealing with beaming angle and bulk Lorentz factor measurements showing a broad range of values.

In this paper we consider the possible ILF forms described above and test these models against observed flux distributions. In order to constrain our model parameters, we use the Fermi GBM and Swift sample of GRBs described in Sect. 2, thereby finding the physical model that explains the data most accurately. In these works, different values of the jet Lorentz factor and beaming angle are considered in order to fit the observed data. In the present paper, we want to test a model where the jets have the same lorentz factor and opening angle and the difference in the observation is due to the offset of the observers and/or to the absorption material along the line of sight. We assume that some fraction of GRBs are observed head-on by observers within the opening angle of the jet (θ ≤ θ₀), while another set of offset observers view the jet from outside the jet opening angle (θ ≥ θ₀). In the second set, the offset observers, we can distinguish between observers just beyond the edge of the jet (θ − θ₀ ≪ θ₀) and those far off-axis (θ − θ₀ ≫ θ₀). Assuming there is no selection bias in favor of larger luminosities, the idea behind the former assumption is that half of all GRBs within a given data set should lie between approximately the 25th and 75th percentile of the luminosity distribution at the redshift of the GRB, and only 25% and 10% should lie below the 25th and the 10th percentile respectively. Therefore, for example, if a sample of four GRBs all lie below the L_median predicted by a particular hypothesis, we can rule out the hypothesis with a confidence of 15/16 or ∼94. If we can compute the value of L_median corresponding to the 25th percentile predicted by a given hypothesis, and all four GRBs lie below this value, then we may rule out the hypothesis to a confidence of (4⁴ − 1)/4⁴ or ∼99.6 %, and so forth. In general, the probability that in a sample of N GRBs all will lie below the jth percentile is (j/100)^N(1 − j/100)^N. Therefore, measurement of as few as for example four GRBs at a given redshift would rule out, with reasonable confidence, a median luminosity of more than ten times their observed average luminosity.

In order to test the model, we first consider a case where the prompt GRB emission is in a narrow(‘pencil’) beam, where all the GRBs are identical, and where the distribution of observers is isotropic. They would nevertheless be distributed in apparent luminosities L because different observers would see them at different viewing angles. If the beam can be considered a pencil beam, then the received luminosity is proportional to 1/(Γ(1 − βcosθ))⁴, whereas the solid angle within the cone defined by the viewing angle is 1 − cosθ. However, for viewing angles that are much greater than 1/Γ, 1 − cosθ ≃ 1 − βcosθ ∝ L^1/4. Hence the value of α should be 5/4.

For GRBs observed slightly beyond the edge of the jet, the source should be considered extended. The effective solid angle from the most visible part of the source is then proportional to the viewing angle from the closest part of the source to the observer’s line of sight, and one picks up an extra factor, namely the Doppler factor. The luminosity therefore roughly follows (1 − cosθ)³, and α takes a value of 4/3. A description of the model is provided by Banerjee et al. (2021). Here, we briefly describe our basic assumptions. We distinguish three possible physical scenarios.

CASE I. Optically thick scenario. The head-on observers do not see any emission as the material in the jet traps the photons, forcing them to exit the jet region, and therefore the jet is optically thick. Observers outside the jet will see different emission according to their position relative to the jet opening angle. As discussed above for viewing angles that are much greater than 1/Γ, 1 − cosθ ≃ 1 − βcosθ ∝ L^1/4. Hence, the value of α should be 5/4. For GRBs observed slightly beyond the edge of the jet, the source should be considered extended. The effective solid angle from the most visible part of the source is then proportional to the viewing angle from the closest part of the source to the observer’s line of sight, and one picks up an extra factor, namely the Doppler factor. The luminosity therefore roughly follows (1 − cosθ)³, and α takes a value of 4/3.

CASE II. Optically thin scenario. For head-on observers, the luminosity function for a set of identical GRBs would be approximately a delta function δ(L − L_max) if the beam were optically thin. This case is excluded from observations because it predicts more bright bursts than are observed. For observers outside the jet, we get the same as CASE I. It is important to notice that this case is excluded because we assume that the prompt emission is produced only within the jet opening angle. However, it may be possible that both short and long GRB jets are structured. In this case, given the jet structure, one can explain the apparent luminosity function of the GRBs by on- and off-axis observers for the optically thin scenario. The low luminosity of GRB 170817A can be explained very well in this model as a prompt emission production above the jet opening angle. It may be that this GRB was detected off axis (see Granot et al. 2018). The luminosity function in the case of the structured jet model was derived by Pescalli et al. (2015; see Eq. (15)). The luminosity function in this case can be represented as a power law with a slope of ∼1.25 (Pescalli et al. 2015). We consider this luminosity function in our analysis (see Table 1).

CASE III. Partially optically thick scenario: The head-on observers may see emission from an optically thin emission region, where the photons last interacted with a region that is obscured from the observer by high optical depth (Eichler 2014). When the shell is optically thick, observers outside the jet see the same LF as the observers outside an optically thin shell because they are seeing photons that are scattered backwards in the frame of the shell. However, observers within the opening angle of the jet cannot see photons backscattered from material moving directly at them. Therefore, there is a ‘blind spot’. These observers see backscattered photons only from jet material outside the blind spot. The calculation of the size of the blind spot, θ_BS, was outlined by Eichler (2014). This latter authors shows that the minimum angle of the blind spot depends on the optical depth of the shell and on the distribution of circumburst material. The blind spot means that there is a minimum offset between the observed material and the line of sight. This is in contrast to observers outside the jet, who can come arbitrarily close to the edge of the jet. Observers within the jet see a light echo from those parts of the jet that are outside the blind spot. Most of the time-integrated signal comes from an annular just outside the blind spot and comparable in solid angle to that of the blind spot itself. If the observer is well within the jet, the annulus completely surrounds the blind spot, but if the blind spot touches the edge of the jet then the annulus is only partial and the observer very close to the edge of the jet sees an annulus that is only about half of the circle. The geometry of this model is discussed in detail by (Vyas et al. 2021). For observers just outside the jet, the luminosity function can be represented by a power law as described above. For observers inside the jet, the luminosity function depends on the distribution of the sizes of blind spots. It is reasonable to represent this latter distribution by another power law. The luminosity function depends on the circumburst material and the illumination of the shell from behind. The total luminosity function is then the sum of the luminosity function for observers outside the jet and the ones inside the jet. Therefore, we always write the ILF as the sum of two terms, N(L)=N_in(L)+N_out(L). The N_out(L) is a power law as discussed above. If the distribution of the blind spot sizes is scale free then it is reasonable to represent N(L)_in as a different power law. In this case, the total ILF can be represented by the sum of two broken power laws: N(L)∝(L/L₀)^−α + ϵN(L)∝(L/L₀)^−β, which can be also approximated as a broken power law, N(L)∝(L/L₀)^−α for L <  L₀ and N(L)∝(L/L₀)^−β for L >  L₀.

3. Median luminosity for gamma-ray bursts

In this section we determine the median luminosity for a given luminosity distribution as a function of the cosmological redshift. In order to evaluate the median luminosity, we use the relation

$$\begin{array}{c}\hfill {\displaystyle {\int }_{L_{\min }}^{L_{\mathrm{median}}}}N(L)\mathrm{d}L={\displaystyle {\int }_{L_{\mathrm{median}}}^{L_{\max }}}N(L)\mathrm{d}L,\end{array}$$(1)

such that L_median is the luminosity dividing the curve into two equal parts. Here, N(L, z) is the luminosity distribution function of a class of emitters, or in other words is the rate of detected GRBs per unit comoving volume at a redshift z between L_min and L_max.

We consider two possible forms of the ILF: The first is the power-law ILF form, N(L)dL ∝ L^−α, which represents the physical case where all the observes are outside the jet θ >  θ₀. We then calculate the median luminosity L_median(z, α) as a function of redshift z.

We show that the data for nearby, low-luminosity GRBs and distant, high-luminosity GRBs are best fitted by a value of α close to 4/3, and that L_median(z) increases with distance by many orders of magnitude within the measured range. Thus, the low luminosity of nearby GRBs is entirely consistent with the universality of GRB energies, and most of the variation in E_iso is attributable to the variation of viewing angles among the observers, who, according to hypothesis, see an ILF index α ranging from 5/4 to 4/3. We note that this hypothesis (Eichler & Levinson 2004, 2006; Levinson & Eichler 2005) neatly explains both the Amati relation (Amati et al. 2002) and the Ghirlanda relation (Ghirlanda et al. 2004).

We consider also a broken power law ILF that represents both the observers inside and outside the jet as explained in the previous section:

$$\begin{array}{c}\hfill N(L)\propto {(L/L_0)}^{-\alpha ,-\beta }.\end{array}$$(2)

Several authors have checked this ILF form and find the best fit values for the free parameters by comparing the expected log N−log S with the observed one (Pescalli et al. 2016; Guetta & Piran 2007). Most of the authors agree that the parameters should be in the range of α ∼ 1.3 (for L <  L₀) and β ∼ 1.8 (for L >  L₀) and a break luminosity L₀ = 10⁵¹ erg s. We tested other possible values of the free parameters using the median luminosity method.

The local luminosity function of GRB peak luminosities L, defined as the comoving space density of GRBs in the interval L to L + dL, is

$$\begin{array}{c}\hfill \mathcal{N}(L)=c_0(\frac{L}{L_{\max }}{)}^{-\alpha },\end{array}$$(3)

where c₀ is a normalization constant with a dimension of 1/L, meaning that the integral over the luminosity function is equal to unity.

In our analysis, we consider the evolution of the luminosity function with redshift z including the k-correction effects, as were considered in our earlier paper (Banerjee et al. 2021). We consider N(L, z)=N(L)N(z) (as in Lloyd-Ronning et al. 2002; Yonetoku et al. 2004), with N(z) defined as (Lloyd-Ronning et al. 2002; Petrosian et al. 2015)

$$\begin{array}{c}\hfill N(z)={\rho }_0\rho (z),\end{array}$$(4)

where the star formation rate is given by Madau & Dickinson (2014)

$$\begin{array}{c}\hfill \rho (z)=\frac{{(1+z)}^{2.7}}{1+{[(1+z)/2.9]}^{5.6}}·\end{array}$$(5)

Here, ρ₀ is the GRB rate at z = 0 expressed in units of Gpc⁻³ yr⁻¹.

Therefore, the evolution of the luminosity function with redshift is given by

$$\begin{array}{c}\hfill N(L,z)={\rho }_0{c}_0(\frac{L}{L_{\max }}{)}^{-\alpha }\frac{{(1+z)}^{2.7}}{1+{[(1+z)/2.9]}^{5.6}}\end{array}$$(6)

out to some maximum redshift z_max, beyond which N is assumed to vanish.

However, for the main focus of this paper, namely the dependence of mean luminosity on redshift, these purely redshift dependence aspects are irrelevant and can be omitted. This can be seen if we substitute the form of N(L, z) from Eq. (6) into Eq. (1). The purely z-dependent terms get cancelled out from right hand side (r.h.s) and left hand side (l.h.s), leaving behind the z dependence coming from L_min as shown below.

Now using Eq. (6) in Eq. (1) and solving for L_median we get,

$$\begin{array}{c}\hfill L_{\mathrm{median}}(z)={\left(\frac{L_{\max }^{-\alpha +1}+{(L_{\min })}^{-\alpha +1}}{2}\right)}^{\frac{1}{-\alpha +1}}.\end{array}$$(7)

Similarly,

$$\begin{array}{c}\hfill L_{25}^{-\alpha +1}-L_{\min }^{-\alpha +1}=\frac{(L_{\max }^{-\alpha +1}-L_{25}^{-\alpha +1})}{3},\end{array}$$(8)

and so

$$\begin{array}{c}\hfill L_{25}=[(\frac{3}{4})(L_{\max }^{-\alpha +1}+L_{\min }^{-\alpha +1}){]}^{1/(1-\alpha )},\end{array}$$(9)

and the luminosity that established the jth percentile is given by

$$\begin{array}{c}\hfill L_j=\frac{L_{\max }^{-\alpha +1}+L_{\min }^{-\alpha +1}}{(1-j)(1-\alpha )}.\end{array}$$(10)

In order to generalize the above expression for any given cosmological model, let us now consider that L_min is established by the detector flux, whereby

$$\begin{array}{c}\hfill L_{\min }=4\pi F_{\mathrm{th}}{D}_{\mathrm{L}}^2.\end{array}$$(11)

Here, D_L is the luminosity distance, which is a function of the redshift z, and F_th is the threshold flux. Therefore, substituting L_min into Eq. (7), we get

$$\begin{array}{c}\hfill L_{\mathrm{median}}(z)={\left(\frac{L_{\max }^{-\alpha +1}+{(4\pi D_{\mathrm{L}}^2(z)F_{\mathrm{th}})}^{-\alpha +1}}{2}\right)}^{\frac{1}{-\alpha +1}}.\end{array}$$(12)

In order to obtain the nature of the variation of L_median, let us consider the Λcold dark matter (CDM) model. The expression for the luminosity distance in terms of the redshift z for the ΛCDM model is given by (assuming the value of the scale factor at the present time is normalized to 1)

$$\begin{array}{c}\hfill D_{\mathrm{L}}=\frac{c}{H_0}(1+z){\displaystyle {\int }_0^{z^{\prime }}}\frac{1}{{({(1+z^{\prime })}^3{\mathrm{\Omega }}_{m_0}+{\mathrm{\Omega }}_{{\mathrm{\Lambda }}_0})}^{1/2}}\mathrm{d}{z}^{\prime },\end{array}$$(13)

where Ω_m0,  Ω_Λ0, and H₀ are the present values of the matter density, the cosmological constant, and the Hubble parameter, respectively.

3.1. Comparing with observed data

We now use Eq. (13) in Eq. (12) and get the variation of L_median with z. However, since the integral can be solved only numerically for ΛCDM, we can only get a numerical solution for the median luminosity.

Observed Data Sample. Most GRBs do not have known redshifts, and knowledge of the redshift, which generally requires identification of a host galaxy, is probably biased towards the high-luminosity GRBs, which receive more attention, and for which it is easier to find a host galaxy (the automatic slewing capabilities of Swift lessen this bias, but do not eliminate it). However, for very nearby GRBs (z ≲ 0.15), a supernova is generally observed, and so no matter how dim the GRB is, its coincidence in time and space with the supernova ensures the determination of its redshift. If there is no associated supernova, coincidence with a nearby Galaxy (which is more statistically significant than coincidence with a distant Galaxy), can reveal the redshift. The discovery of a kilonova or a gravitational wave signal can, as in the case of GRB 170817 (LIGO Scientific Collaboration & Virgo Collaboration 2017a,b,c; Goldstein et al. 2017; Savchenko et al. 2017) help in confirming the host galaxy. LIGO and Virgo observations can shed light on short GRBs, but not the mostly long GRBs that we have considered for our analysis. In our present work, the observed data are taken from Goldstein et al. (2017), Cano et al. (2017), Michałowski et al. (2018), Gehrels et al. (2006), D’Avanzo et al. (2014), Demianski et al. (2017). Our sample consists of data collected from different sources in such a way that it spans a wide range of redshift starting from the nearest GRB (z = 0.00867) up to redshifts as high as 9. The low-redhsift GRBs are mostly associated with supernovae. The sample has also been chosen such that we have adequate data (including the values of spectral peak, redshift, and flux) for each GRB so that they can be included in all our plots. The Swift data comprise 17 short GRBs, 115 long GRBs, and 3 ultra-long GRBs.

We further divided the data sets based on their afterglow nature: flat or monotonic.

In Banerjee et al. (2021), we obtained the best-fit parameter values for single and broken power-law luminosity functions using peak flux distribution data for Swift and Fermi GRBs. The Fermi GRB sample consisted of 918 long GRBs obtained from Gruber et al. (2014), Kienlin et al. (2014), Bhat et al. (2016), Meszaros & Meszaros (1996). The Swift sample on the other hand consisted of 365 long GRBs obtained from the Swift catalogue, Gruber (2011, 2012). For both cases, we obtained the histogram for the differential source counts using their peak flux distributions and compared it to the peak flux of the whole sample with known and unknown redshifts. Using the method of χ² minimization and p value, we obtained best-fit parameter values of (i) α = 1.17 for single power law and (ii) $\alpha =0.6_{-1}^{+1},\beta =1.{25}_{-0.05}^{+0.05}$, and $L_0=5_{-2}^{+2}\times {10}^{51}$ for a broken power law.

In Fig. 1, we plotted the variation of L_median with redshift for α = 1.5,  1.33,  1.25, 1.41, and 1.17, and L_max = 3 × 10⁵³ for two different values of threshold flux F_th = 2 × 10⁻⁸ and F_th = 1 × 10⁻⁷. We also show the results for L_median corresponding to the ILF given by a broken power law defined in Eq. (2) using the best-fit parameter values obtained from our analysis of the log N − log F plot in Banerjee et al. (2021). Here, we show that the median luminosity can also be used to check whether or not these models are accurate. Henceforth, L and F correspond to the peak bolometric luminosity and flux, respectively, with L in units of erg s⁻¹ and F_th in units of erg cm⁻² s⁻¹.

Fig. 1. Plot showing the variation of Lmedian (in erg s−1) with z along with observed Swift data and simulated data. The plots in panel a correspond to Fth = 2 × 10−8 erg cm−2 s−1. The plots in panel b correspond to Fth = 1 × 10−7 erg cm−2 s−1. The values of the parameters are marked on the plots. The orange dots represent the simulated data. The other, larger dots correspond to the observed data for long, short, and ultralong GRBs (Goldstein et al. 2017; Cano et al. 2017; Michałowski et al. 2018; Gehrels et al. 2006; D’Avanzo et al. 2014; Demianski et al. 2017). The terms monotonic, flat, and unavailable correspond to the nature of afterglow data for the corresponding GRB. The curve marked threshold represents the threshold value (for the rest of the plots in this paper) of L at a given z. For comparison, we also plot the curve (green and brown) corresponding to the ILF given by a broken power law.

In Fig. 2, we highlight the low-z region. We also show the position of GRB170817 (Goldstein et al. 2017) with L_peak = 1.3 × 10⁴⁷erg s⁻¹ and z = 0.009 in our plot. In both figures, the values of α are mentioned in the plot as well.

Fig. 2. Variation of Lmedian (in erg s−1) with z for small redshifts and comparison with simulated and observed data. Fth is in units of erg cm−2 s−1. The plots in the panel a correspond to Fth = 2 × 10−8. The plots in panel b correspond to Fth = 1 × 10−7. The values of the spectral index α are marked on each curve. The orange dots are for simulated data. The larger dots correspond to the observed data for long, short, and ultralong GRBs (Goldstein et al. 2017; Cano et al. 2017; Michałowski et al. 2018; Gehrels et al. 2006; D’Avanzo et al. 2014; Demianski et al. 2017). The green dots are for short GRBs with monotonic afterglow, magenta dots for long GRBs with unavailable afterglow data, and the yellow dots correspond to ultralong GRBs with monotonic afterglow data.

As we are plotting L_median, the curve that best fits the data should have an equal number of data points below and above the corresponding theoretical curve. Counting the number of data points above and below each curve for a single power law in Fig. 1 we find that for F_th = 2 × 10⁻⁸, the curve corresponding to α = 1.17 is the most symmetrical one, with α = 1.25 lying close by. The curves corresponding to other values of α are much less symmetrical and can therefore be rejected. We find that if we increase the value of F_th by a factor of five as shown in Fig. 1b, the most symmetrically curve for a single power law is now given by a higher value of α = 1.41 with α = 1.33 lying close by.

For the broken power law, we plot the best-fit values obtained in Banerjee et al. (2021): α = 0.6,  β = 1.25, and L₀ = 5 × 10⁵¹. We find that a broken power law with these parameter values fits the data most closely for F_th = 2 × 10⁻⁸. For F_th = 1 × 10⁻⁷, we obtain the same results. The exact numbers of the counts are shown in Tables 1 and 2.

For very low redshifts, as in Fig. 2, the broken power law reduces to a single power law case, depending only on the spectral index α. Hence, we do not show it in Fig. 2 as it is not relevant.

For the present case, in order to quantify the goodness of the fit using the present sample, we compare the theoretical and observed results to obtain the p values for single and broken power laws using the best-fit parameter values obtained in Banerjee et al. (2021). For this, we divided the analysis into ten redshift bins, and compared the median luminosity from the observed data and the theoretical median luminosity for each redshift bin. In Table 3, we show the p values obtained using the best-fit parameter values for both single and broken power law cases.

We note that we also plot the threshold curve in all the figures for reference. Some of the data points lie either on or slightly below it. However, these GRBs are either very long, have uncertainties in their redshifts, or are affected by X-ray flash. Hence, all of them have properties different from other GRBs. From Fig. 1b, which contains the Swift data, we see that, among the long GRBs, the ones with flat afterglows are much less scattered compared to those with monotonic afterglows. The probability of this being by chance is 0.03 or 3% at 95% confidence. Thus, this difference is statistically significant.

We further checked whether or not lowering L_max causes any improvement. We find that the change in L_median due to a change in L_max by a factor of ten is almost negligible. Hence, introducing a spread in L_max given by a sum of two luminosity functions will not cause any significant change in L_median.

As we know, for head-on observers, the ILF for a set of identical GRBs is a delta function δ(L − L_max) if the beam is optically thin. Following our earlier analysis, one can easily see that the L_median for this case is a constant around L_max which clearly does not match the observational data, and is therefore ruled out.

3.2. Comparing with simulated data

Observed data may be biased, especially at high redshifts. Hence, in order to improve our statistics, we further compared our theoretical results with simulated data for L_median vs. z. In order to generate the simulated data, we considered a pencil beam with random orientations of GRBs in the sky. In other words, we assumed that GRBs can occur anywhere within a certain volume at equal cosmic times. We further divided intervals of redshift into equal volumes. As discussed above, the dependence of luminosity on the jet opening angle for a pencil beam is proportional to 1/(Γ(1 − βcosθ)⁴, θ being the jet opening angle.

Using this information, we generated the simulated data for L_median in terms of redshift with an average of 17 GRBs in each redshift bin(dz = 0.01). The results are shown in Figs. 1 and 2 for both large and small redshifts for the two threshold values we considered in the earlier section.

In order to compare the theoretical results with simulated data, we once again counted the number of points above and below each theoretical curve and identified the curve that divides the data set most symmetrically. We repeat this for both Figs. 1a and b. As above, we find that for F_th = 2 × 10⁻⁸, among the chosen values of α = 1.17,  1.25,  1.33,  1.41,  and 1.5, the L_median curves for α = 1.17 divide the data points most symmetrically with α = 1.25 lying close by. Higher values of α worsens the fit. As we increase the value of F_th to 1 × 10⁻⁷, the best-fit curve for a single power law is now given by α = 1.33. The exact numbers are shown in Tables 1 and 2.

For completeness, using the best-fit values of the free parameters obtained above, we see that they once again divide the simulated points most symmetrically. The exact numbers, for both the thresholds, are shown in Tables 1 and 2.

Thus, comparing our results with both simulated and observed data, we can safely conclude that for the case where the luminosity of the GRBs is given by a power law only, their luminosity distribution with respect to redshift is best given by a pencil beam, that is, α = 1.17 and α = 1.25. Higher values of α fit the data very poorly, and are therefore rejected. This result is strong evidence that the jet may be structured as shown in Pescalli et al. (2015). The best-fit values of the broken power-law parameters given by Banerjee et al. (2021) give a very good fit to the data for the present case also.

3.3. Analytical solution for an extended jet

Let us now consider an extended jet with an opening angle θ₀ such that dL₀/dΩ′ = δ(θ₀ − θ). Here, dL₀ is the contribution from the head-on element dΩ′ of the jet, and θ is the angle between the jet axis (situated at zero) and the observer. Considering the Doppler effect,

$$\begin{array}{c}\hfill \frac{\mathrm{d}L}{\mathrm{d}{\mathrm{\Omega }}^{\prime }}=\frac{\mathrm{d}{L}_0}{\mathrm{d}{\mathrm{\Omega }}^{\prime }}\frac{{(1-\beta )}^4}{{(1-\beta cos\chi )}^4},\end{array}$$(14)

where χ is the angle between the observer and the element dΩ′, the total luminosity is obtained by integrating the above expression over dΩ′:

$$\begin{array}{c}\hfill L={\displaystyle \int \frac{\mathrm{d}L}{\mathrm{d}{\mathrm{\Omega }}^{\prime }}}\mathrm{d}{\mathrm{\Omega }}^{\prime }={\displaystyle \int \frac{\mathrm{d}{L}_0}{\mathrm{d}{\mathrm{\Omega }}^{\prime }}}\frac{{(1-\beta )}^4}{{(1-\beta cos\chi )}^4}\mathrm{d}{\mathrm{\Omega }}^{\prime }.\end{array}$$(15)

Writing cosχ = sinθsinθ′ cos ϕ′+cos θcos θ′, we get

$$\begin{array}{c}\hfill L={\displaystyle \int \delta }({\theta }_0-\theta )\frac{{(1-\beta )}^{4}sin{\theta }^{\prime }\mathrm{d}{\theta }^{\prime }\mathrm{d}{\varphi }^{\prime }}{{(1-\beta (sin\theta sin{\theta }^{\prime }cos{\varphi }^{\prime }+cos\theta cos{\theta }^{\prime }))}^4}.\end{array}$$(16)

The luminosity distribution function is now given by

$$\begin{array}{c}\hfill \mathrm{d}N=\frac{\frac{\mathrm{d}N}{\mathrm{d}cos\theta }}{\frac{\mathrm{d}L}{\mathrm{d}cos\theta }}\mathrm{d}L.\end{array}$$(17)

As we know, dN/dcos θ = constant = K where K is the rate of GRBs, which is two per day.

Using the above ILF, one can now compute the median luminosity, and therefore find the best-fit values of Γ and θ₀. For this, we use the Fermi GBM data as before. In Fig. 3a, we plotted the variation of L_median with z, setting θ₀ = 0.1. Different curves correspond to different values of Γ. Plotting the curves along with GBM data, we find from Fig. 3a that for θ₀ = 0.1, the best-fit value of Γ = 5.15. Γ = 100 also gives a good fit. However, values of Γ such as 10.01,  31.6, and so on give a very poor fit, and Γ values of < 2.7 fit the data quite well. In Fig. 3b, we plot the curves for different values of θ₀ for Γ = 100. Comparing with data, we find that the best-fit value of θ₀ is 0.1, which is already clear from Fig. 3a. We also see that values of θ₀ of lower than 0.1 still give a good fit, but the fit becomes poorer above θ₀ = 0.1.

Fig. 3. Variation of Lmedian with z for an extended jet along with Fermi GBM data. The plots in panel a are for θ = 0.1. Different curves correspond to different values of β. Curves in panel b are for β = 0.99995.

In Fig. 4, we compare the numerical solution of Eq. (16) to analytical approximate expressions [1 + (θ − θ_o)²(Γ)]⁻³ and [1 + (θ)²(Γ)²]⁻⁴. In Fig. 5, we show a simulated plot for L weighted by a star formation z dependence (1 + z)^k for k = 1.4 (Lloyd-Ronning et al. 2002; Yonetoku et al. 2004). This plot was obtained for θ₀ = 0.1,  Γ = 100, and L_max = 10⁵³. Points are above the Fermi GBM threshold are shown. We also show the threshold curve and analytical curves for α = 1.33 and 1.25. As we can see, the contribution of jet elements outside the opening angle – that is, larger values of θ − θ₀ – becomes important and observable only near low z. For higher redshifts, the main contribution comes from elements near to the jet axis, that is, in and around the opening angle. We also show the real L_median values that divide the number of GRBs equally at a given z. In order to clearly see the effect of including the GRBs within the jet, we show the real L_median for both with and without the GRBs inside the jet for Γ = 100. As we can see, when we include the GRBs within the jet, the L_median curve no longer passes through the data point, and therefore gives a poor fit. However, for highly opaque medium, that is, where we can see only those parts of the jet where the velocity vector is sufficiently misaligned with the line of sight, the L_median curve gives a much better fit. For reference, we also show L_median for Γ = 30,  Γ = 200.

Fig. 4. Numerical solution of Eq. (16) compared to analytical approximate expressions: [1 + (θ − θ0)2(Γ)]−3 (in the region just beyond θ0) represented by the green dotted line and [1 + (θ)2(Γ)2]−4 (the region 0.2 <  θ <  0.3) represented by blue dashed line.

Fig. 5. Simulated data showing variation of L with z for different ranges of θ. The black stars correspond to the value of L(z) that divides the points equally (above and below). The green stars reflect the same L(z) but without the green points, i.e., for points outside the jet. The yellow squares and cyan rectangles represent similar points (L(z)) for Γ = 30 and Γ = 200, respectively. The other, larger diamonds correspond to the observed data (Goldstein et al. 2017; Cano et al. 2017; Michałowski et al. 2018; Gehrels et al. 2006; D’Avanzo et al. 2014) along with the nature of their afterglows (flat, monotonic, or unavailable): dark cyan = short, flat; dark magenta = flat, monotonic; dark yellow = flat, unavailable; dark pink = long, flat; black = long, monotonic; and dark orange = long, unavailable. The analytical curves for α = 1.25 (black dot-dashed) and α = 1.33 (black dashed) are also shown.

Finally, in Fig. 6, we plot the E_p vs. L_iso, varying both with respect to θ with a spread in L_max ∼ 𝒪(10). From Eichler & Levinson (2006) and Ioka & Nakamura (2001), we know that the maximum frequency goes as ν_max ∝ 1/(1 − βcos χ). Therefore, expressing χ in terms of θ and integrating over the jet element, we get E_p as a function of θ. Figure 6 was obtained for θ₀ = 0.1 and Γ = 100. We plot $E_{\mathrm{p}}\propto L_{\mathrm{iso}}^{1/4}$ and $E_{\mathrm{p}}\propto L_{\mathrm{iso}}^{1/3}$ curves for reference. As we can see, for larger angles (away from the edge), the simulated plot follows the curve $E_{\mathrm{p}}\propto L_{\mathrm{iso}}^{1/4}$ and then near the edge it goes more like $E_{\mathrm{p}}\propto L_{\mathrm{iso}}^{1/3}$.

Fig. 6. Ep vs. Liso as a function of θ along with observed Fermi GBM data. The ranges of θ are the same as in Fig. 5. The lines corresponding to $E_{\mathrm{p}}\propto L_{\mathrm{iso}}^{1/4}$ and $E_{\mathrm{p}}\propto L_{\mathrm{iso}}^{1/3}$ are also shown.

4. Conclusions

In this paper, we study the behaviour of the quantity L_median for a given luminosity function. We studied the behaviour for the luminosity functions given by a power law and we find that, under very general assumptions, L_median always increases with redshift and therefore distance, as expected. We also find it to be very insensitive to L_max, the intrinsic luminosity seen by an on-axis observer. In our present work, we have chosen a class of physical models described by Eichler (2014; and references therein). This model is essentially a shell with baryonic material illuminated from behind that is possibly optically thick. For simplicity, we consider a constant bulk Lorentz factor and an opening angle θ₀. The GRB ILF consists of a contribution of N_out(L) observers outside the jet and N_in(L) inside θ₀. We calculated N_out(L) and show that it can be approximated as N(L)dL ∝ L^−αdL where α = 4/3 for observers just outside the jet and α = 5/4 for observers far from the jet. This power law extends to a maximum value that corresponds to the observers within 1/Γ of the edge of the jet. For observers inside the jet, N_in(L) depends on the distribution of blind spot sizes because the optical depth is large when the shell is illuminated from behind. In this case, the blind spot hides material within it. It was argued by Eichler (2014) that the blind spots size is likely to be larger than 1/Γ. Therefore, the brightest bursts may well be dominated by N_out(L) and possibly some optically thin bursts. If N_in(L) is represented by a power law, then it should be softer, more biased towards lower luminosities, and have a steeper power-law index. If at the other extreme the blind spots have all the same size, then N_in(L) is represented by a delta function that peaks at a much lower luminosity than the maximum luminosity. As it is difficult to calculate the distribution of the sizes of blind spots, we considered both extremes. For completeness, we also considered the possibility that the GRB jet is structured and considered the implied luminosity function found by Pescalli et al. (2015). We find that this model can represent the observed data very well. For this, we compared the median luminosities obtained theoretically and from the observed data using the best-fit parameter values obtained by Banerjee et al. (2021). The results are presented in Table 3. As we can see, the fit is good even for the present case with median luminosity. We find that the best value of α calculated in this paper and the best-fit value of β and the luminosity break are in agreement with the results of Banerjee et al. (2021) and with other values found in the literature.

Acknowledgments

We would like to acknowledge the contribution of Prof. David Eichler towards the development of this work, who sadly passed away during the elaboration of this manuscript. We would like to thank with Avishai Gal-Yam, Richard Ellis, Elena Pian and Paolo Mazzali for useful discussions. We acknowledge funding from Israel Science Foundation and John and Robert Arnow Chair of Theoretical Astrophysics.

References

All Tables

All Figures

Fig. 1. Plot showing the variation of Lmedian (in erg s−1) with z along with observed Swift data and simulated data. The plots in panel a correspond to Fth = 2 × 10−8 erg cm−2 s−1. The plots in panel b correspond to Fth = 1 × 10−7 erg cm−2 s−1. The values of the parameters are marked on the plots. The orange dots represent the simulated data. The other, larger dots correspond to the observed data for long, short, and ultralong GRBs (Goldstein et al. 2017; Cano et al. 2017; Michałowski et al. 2018; Gehrels et al. 2006; D’Avanzo et al. 2014; Demianski et al. 2017). The terms monotonic, flat, and unavailable correspond to the nature of afterglow data for the corresponding GRB. The curve marked threshold represents the threshold value (for the rest of the plots in this paper) of L at a given z. For comparison, we also plot the curve (green and brown) corresponding to the ILF given by a broken power law.

In the text

Fig. 2. Variation of Lmedian (in erg s−1) with z for small redshifts and comparison with simulated and observed data. Fth is in units of erg cm−2 s−1. The plots in the panel a correspond to Fth = 2 × 10−8. The plots in panel b correspond to Fth = 1 × 10−7. The values of the spectral index α are marked on each curve. The orange dots are for simulated data. The larger dots correspond to the observed data for long, short, and ultralong GRBs (Goldstein et al. 2017; Cano et al. 2017; Michałowski et al. 2018; Gehrels et al. 2006; D’Avanzo et al. 2014; Demianski et al. 2017). The green dots are for short GRBs with monotonic afterglow, magenta dots for long GRBs with unavailable afterglow data, and the yellow dots correspond to ultralong GRBs with monotonic afterglow data.

In the text

	Fig. 3. Variation of Lmedian with z for an extended jet along with Fermi GBM data. The plots in panel a are for θ = 0.1. Different curves correspond to different values of β. Curves in panel b are for β = 0.99995.
In the text

	Fig. 4. Numerical solution of Eq. (16) compared to analytical approximate expressions: [1 + (θ − θ0)2(Γ)]−3 (in the region just beyond θ0) represented by the green dotted line and [1 + (θ)2(Γ)2]−4 (the region 0.2 <  θ <  0.3) represented by blue dashed line.
In the text

Fig. 5. Simulated data showing variation of L with z for different ranges of θ. The black stars correspond to the value of L(z) that divides the points equally (above and below). The green stars reflect the same L(z) but without the green points, i.e., for points outside the jet. The yellow squares and cyan rectangles represent similar points (L(z)) for Γ = 30 and Γ = 200, respectively. The other, larger diamonds correspond to the observed data (Goldstein et al. 2017; Cano et al. 2017; Michałowski et al. 2018; Gehrels et al. 2006; D’Avanzo et al. 2014) along with the nature of their afterglows (flat, monotonic, or unavailable): dark cyan = short, flat; dark magenta = flat, monotonic; dark yellow = flat, unavailable; dark pink = long, flat; black = long, monotonic; and dark orange = long, unavailable. The analytical curves for α = 1.25 (black dot-dashed) and α = 1.33 (black dashed) are also shown.

In the text

	Fig. 6. Ep vs. Liso as a function of θ along with observed Fermi GBM data. The ranges of θ are the same as in Fig. 5. The lines corresponding to $E_{\mathrm{p}}\propto L_{\mathrm{iso}}^{1/4}$ and $E_{\mathrm{p}}\propto L_{\mathrm{iso}}^{1/3}$ are also shown.
In the text