© ESO, 2013

1. Introduction

Gamma-ray bursts (GRBs) are among the most violent and energetic phenomena observed in the Universe. The total energy emitted (if isotropically) in the form of energetic photons can attain values of up to few times 10⁵⁴ erg, two orders of magnitude higher than the energy released by supernova explosions. However, this dramatic energy output can be considerably alleviated if the emission is strongly collimated along jets (Rhoads 1999).

GRBs are classified into two main classes: long (LGRB) and short (SGRB) (Kouveliotou et al. 1993), according to the observed burst duration. The first category includes events with durations longer than ~2 s, while the second includes those with durations shorter than ~2 s. More recent investigations indicate that there may be an intermediate class of bursts that has the softest energy spectra of the three (Horváth et al. 2008, 2010; Huja et al. 2009). The comparison between Swift and BATSE samples indicates that the population fraction in these classes differs because Swift detectors are more sensitive to soft and weak bursts. As a consequence, the probability detection of bursts included in the intermediate group is enhanced compared with the BATSE results (Horváth et al. 2010). Considering only the more conservative classification into two categories, Swift data indicate that 83% of the events are LGRBs, while only 17% are included in the SGRB class (Huja et al. 2009). It is worth mentioning that these two classes have indeed some distinct properties, which are examined below.

For instance, in the plane defined by the temporal lag between features in the light curve observed at different wavelengths and the peak luminosity, these two aforementioned classes are clearly distributed separately. Moreover, there are other indications that LGRBs and SGRBs represent events with different origins. LGRBs are generally located in active star formation regions preferentially found in low-metallicity and low-luminosity galaxies (Fruchter et al. 2006). These events might be related with type Ibc supernovae because several unambiguous spectroscopic identifications were obtained up to now (Stanek et al. 2003; Hjorth et al. 2003; Woosley & Bloom 2006; Berger et al. 2011). However, late radio observations of 68 SNIbc, searching for emission attributable to off-axis events, indicate that less than 10% of type Ibc supernovae might be related to LGRBs (Soderberg et al. 2006b). Recent Swift data indicate that SGRBs occur in galaxies with different morphologies, including early-type objects in which no star formation activity is observed (Berger 2011; Church et al. 2011). Moreover the positioning of these bursts indicates a wide range of spatial offsets from the (assumed) host galaxy (Church et al. 2011). This would be expected if SGRBs were the consequence of a merger of two compact objects that have received a kick velocity during their formation (Church et al. 2011). Moreover, as mentioned above, in addition to the possible existence of an intermediate class, some authors have suggested two other subclasses: subluminous long GRBs (SL-LGRBs) and short GRBs with extended emission (SGRB-EE; Coward 2005; Murase et al. 2006; Guetta & Della Valle 2007; Imerito et al. 2008). But, possible criticisms not withstanding these subclasses are ignored in our investigation and only LGRBs are considered. If LGRBs are the consequence of the death of massive stars, their formation rate R_grb as a function of the redshift is expected to be proportional to the cosmic star formation rate (CSFR) R_∗(z). The observed rate of LGRBs, which permits estimating the local formation rate, will depend not only on the CSFR, but also on the luminosity function φ(L). According to Salvaterra & Chincarini (2007), models in which LGRBs trace the CSFR and are described by a non-evolving luminosity function underestimate the number of events at high redshift. This difficulty could be surpassed if the luminosity function evolves in the sense that LGRBs are more luminous at high redshift (Butler et al. 2010). Including a moderate evolution in the luminosity function and considering an isotropic emission model, Salvaterra & Chincarini (2007) derived a local LGRB formation rate R_grb = 0.12 Gpc^-3 yr^-1, while more recently, Salvaterra et al. (2011) argued that a strong evolution of the luminosity function is required to explain the redshift distribution of these events. A different method was adopted by Wanderman & Piran (2010), who derived the luminosity function and the local formation rate by directly inverting the redshift distribution of LGRBs. They estimated a local rate R_grb = 1.3 Gpc^-3 yr^-1, increasing up to z ~ 3 and decreasing for high redshift. Another possibility that could explain the excess of events at high redshift is to assume that the efficiency of massive stars of forming GRBs increases with redshift. A variable efficiency of forming GRBs could be an effect due to the metallicity of the environment (Daigne et al. 2006; Salvaterra & Chincarini 2007; Qin et al. 2010; Cao et al. 2011; Virgili et al. 2011). If one accepts that LGRBs are formed preferentially in a metal-poor medium, the association of (at least a fraction of them) with type Ibc supernovae is difficult to understand since SNIbc preferentially occur in galaxies with high metallicity (Prantzos & Boissier 2003; Prieto et al. 2008). Pélangeon et al. (2008) followed a different approach to estimate the local formation rate of GRBs. They calculated the volume detectability for sources listed in the HETE-2 catalog, that is the volume V_max of the Universe in which events are bright enough to be included in their sample. Then, the number of GRBs inside V_max was calculated, which permits an estimate of R_grb. They obtained, before correcting for beaming, a local formation rate R_grb ~ 11 Gpc^-3 yr^-1, a value substantially higher than the estimates mentioned previously.

In the present work a novel approach was adopted to derive the local formation rate of LGRBs and the redshift distribution of these events. A Monte Carlo code was developed to estimate the probability that a GRB generated at a given redshift is detected in a given energy band of the detector (visibility function) which in our case corresponds to the Swift-BAT experiment. Once the visibility function is known, the formation rate and the redshift distribution of GRBs can be estimated. Another particular aspect of our approach is that instead of introducing the luminosity function of GRBs as is commonly done, we consider the energy distribution of the events, that is the probability for a given burst to occur with an energy E. The best parameters defining the assumed energy distribution were then determined by comparing predicted fluences with those given by the Swift-BAT catalog. The choice of fluences is dictated by the fact that this primary physical quantity is well defined and well measured. The paper is organized as follows: in Sect. 2 the jet model and the Monte Carlo code are described; in Sect. 3 we discuss the simulations and the data analysis that lead to an estimate of the parameters that define the energy distribution; in Sect. 4 the main results are given, and finally we present the conclusions in Sect. 5.

2. Model and the Monte Carlo code

Two phenomenological jet models have been intensively discussed in the literature. The so-called universal model, in which jets have the same structure but the angular energy distribution varies (Meszaros et al. 1999; Rossi et al. 2002) and the alternative model, which considers that the angular energy distribution is uniform inside the jet, but the total emitted energy is nearly the same for all events, although beamed into different opening angles (Rhoads 1997; Dai & Cheng 2001; Frail et al. 2001).

In our approach, we assumed that the γ-ray emission occurs along two cones of aperture 2θ. In this case, the total energy E_J emitted along the jets is related to the isotropic energy E_iso by the well-known relation ${\mathit{E}}_{\mathrm{J}}\mathrm{=}{\mathit{E}}_{\mathrm{iso}}\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathit{\theta }\mathrm{)}\mathit{.}$(1)We assumed furthermore that the probability for an event to occur with an energy E_J is given by a log-normal distribution, that is $\mathit{P}\mathrm{\left(}{\mathit{E}}_{\mathrm{J}}\mathrm{\right)}\mathrm{d}\ln {\mathit{E}}_{\mathrm{J}}\mathrm{=}\mathit{A}\exp \left[\mathrm{-}\frac{{\left(\ln {\mathit{E}}_{\mathrm{J}}\mathrm{-}\ln {\mathit{E}}_{\mathrm{med}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\ln {\mathit{E}}_{\mathrm{j}}}^{\mathrm{2}}}\right]\mathrm{d}\ln {\mathit{E}}_{\mathrm{J}}\mathit{,}$(2)where A is a normalization constant. The median, lnE_med, and the dispersion, σ_lnEJ, of the log-normal distribution are free parameters to be determined. These parameters are constants if there is no evolution in the energy distribution. Since this possibility cannot be excluded, cases in which the median energy evolves were also considered and were modeled simply by the relation $\ln {\mathit{E}}_{\mathrm{med}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}\ln {\mathit{E}}_{\mathit{o}}\mathrm{+}\mathit{\delta }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}\mathit{,}$(3)where lnE_o and δ are constants (note that the median energy of the distribution at z = 0 is E_med(0) = E_oe^δ). The assumption that at high redshift LGRBs might be more energetic may be justified by the following argument. If jets are related to the rotation of the GRB progenitor, it is expected that in the past, massive stars were formed preferentially in a metal-poor environment and, consequently, these stars have mass-loss rates lower than stars formed in metal-rich regions. Less massive axisymmetric winds carry out less angular momentum (de Araujo & de Freitas Pacheco 1994; de Araujo et al. 1994). In this case, the progenitors of LGRBs formed in a metal-poor environment are expected to have higher rotation rates and to produce stronger jets. Note that in this context the metallicity is expected to affect the energetics of the event, but not necessarily the formation efficiency. We cannot exclude the possibility that both the efficiency to produce GRBs and the energy (or luminosity) distribution vary with redshift, but in this investigation we considered a constant efficiency per unit of stellar mass formed and possible evolutionary effects only in the energy distribution, as modeled by Eq. (3).

We assumed also that the probability G(θ) for having a jet with an opening half-angle θ is independent of the jet energy. In other words, the probabilities P(E_J) and G(θ) are statistically independent. The observed distribution of opening angles Q(θ) is related to the intrinsic distribution by the relation $\mathit{Q}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{d}\mathit{\theta }\mathrm{=}\mathit{G}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{(}\mathrm{1}\mathrm{-}\cos {\mathit{\theta }}^{\mathrm{)}}\mathrm{d}\mathit{\theta }\mathrm{\simeq }\frac{\mathrm{1}}{\mathrm{2}}{\mathit{\theta }}^{\mathrm{2}}\mathit{G}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{d}\mathit{\theta }$(4)where we used the fact that the probability for the jet with an opening half-angle θ to be aligned with the line of sight is (1−cosθ). From Eq. (4) the intrinsic distribution G(θ) can be derived if the observed distribution Q(θ) is known.

Estimates of the opening angle of the jet depend on the model parameters. Racusin et al. (2009) computed opening angles for 28 GRBS based on the analysis of the break time of the conical blast wave describing X-ray afterglows. A large sample was considered by Goldstein et al. (2011), who used the Ghirlanda lower limit (Ghirlanda et al. 2004) in the plane ε_peak − fluence, calibrated with GRBs that have well-constrained jet opening angles. These studies suggest that the observed distribution of values of θ can be represented by a log-normal distribution. Although this may be only a consequence of the assumptions made for estimating θ, we adopted this in our computations, that is $\mathit{Q}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{d}\ln \mathit{\theta }\mathrm{=}\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }{\mathit{\sigma }}^{\mathrm{2}}}}\exp \left[\mathrm{-}\frac{{\left(\ln \mathit{\theta }\mathrm{-}\ln {\mathit{\theta }}_{\mathrm{med}}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]\mathrm{d}\ln \mathit{\theta }\mathit{.}$(5)Our preferred model has an observed median value for the jet opening angle 2θ_med = 16° and a dispersion σ = 0.76. It should be mentioned that recently Lu et al. (2012) claimed that the opening angles are anti-correlated with the redshift, meaning that jets were narrower at high z. This possibility is not considered here although it deserves more detailed investigation.

2.1. Code

Our code is similar to that developed by Lamb et al. (2005), who aimed to investigate the unified jet model. For the sake of completeness, we present here the main steps defining the sequence of calculations that permitted us to produce a mock catalog of LGRBs, from which one can identify the differences between the present approach and that of the aforementioned authors.

• 1)

  Initially, the redshift of the object is determined by a probability function proportional to the CSFR. The adopted CSFR was parametrized by the usual form proposed by Cole et al. (2001), that is ${\mathit{R}}_{\mathrm{\ast }}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathit{h}}_{\mathrm{70}}\frac{\mathrm{(}\mathit{a}\mathrm{+}\mathit{bz}\mathrm{)}}{\left[\mathrm{1}\mathrm{+}\mathrm{(}\mathit{z}\mathit{/}\mathit{c}{\mathrm{)}}^{\mathrm{d}}\right]}\mathit{,}$(6)where a = 0.0103, b = 0.12, c = 5.0, d = 2.8, with the CSFR given in M_⊙ Mpc^-3 yr^-1. h₇₀ is the Hubble parameter in units of 70 km s^-1 Mpc^-1. This expression provides a quite good fit of the existing data (see, for instance, Hopkins & Beacom 2006) and is consistent with the cosmological simulations by Filloux et al. (2010, 2011). Under these conditions, the probability Ψ(z) of generating a burst in the redshift interval z,z + dz is $\mathrm{\Psi }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{d}\mathit{z}\mathrm{=}{\mathit{N}}_{\mathrm{\ast }}^{-1}\frac{{\mathit{R}}_{\mathrm{\ast }}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\frac{\mathrm{d}\mathit{V}}{\mathrm{d}\mathit{z}}\mathrm{d}\mathit{z,}$(7)where dV is the comoving volume differential element and the normalization constant is defined by ${\mathit{N}}_{\mathrm{\ast }}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{z}}_{\max }}\frac{{\mathit{R}}_{\mathrm{\ast }}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\frac{\mathrm{d}\mathit{V}}{\mathrm{d}\mathit{z}}\mathrm{d}\mathit{z,}$(8)where we assumed z_max = 15, but the results are not significantly modified for higher values. After obtaining the redshift z from the probability function Ψ(z), the luminosity distance D_L is computed from the adopted cosmology which, in our case, was a ΛCDM model defined by the parameters H₀ = 70 km s^-1 Mpc^-1, Ω_m = 0.3 and Ω_v = 0.7.

• 2)

  In a second step, the energy of the jet is obtained by using the probability function P(E_J) for given values of the median and the dispersion, according to cases including or excluding the evolution of the burst mean energy.

• 3)

  The opening angle of the jet was derived from the probability function G(θ) as described previously.

• 4)

  In the next step, we determine whether the jet points in the direction of the observer with a probability (1−cosθ). If the jet is not visible, the event is saved with its corresponding redshift as an unseen burst. If the event occurs in the observer’s direction, we still have to verify whether the resulting photon rate in the observer’s frame will be able to trigger the considered detector.

• 5)

  We assumed that the photon energy distribution of the burst can be represented by a Band function (Band et al. 1993) with exponents α = −1.0 and β = −2.25, typical values found from fits of existing data. In this case, for photon energies ε ≤ 1.25ε_peak the photon spectrum is $\mathit{B}\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}\mathrm{=}\mathit{K}{\mathit{\epsilon }}^{-1}\exp \mathrm{(}\mathrm{-}\mathit{\epsilon }\mathit{/}{\mathit{\epsilon }}_{\mathrm{peak}}\mathrm{)}\mathit{,}$(9)where K is a normalization constant. For photon energies ε > 1.25ε_peak the photon energy distribution is $\mathit{B}\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}\mathrm{=}\mathrm{0.3786}\mathit{K}{\mathit{\epsilon }}_{\mathrm{peak}}^{\mathrm{1.25}}{\mathit{\epsilon }}^{-2.25}\mathit{.}$(10)In these relations ε_peak is the energy at which the function ε²B(ε) attains a maximum. The value of ε_peak for a given simulated event was estimated from the so-called Amati relation, which correlates ε_peak with the equivalent isotropic energy of the burst. In our simulations E_iso was estimated from Eq. (1), using E_J obtained in step (2). From data by Amati et al. (2008, 2009) the following fit was obtained and adopted for our simulations (rest frame quantities): $\log {\mathit{\epsilon }}_{\mathrm{peak}}\mathrm{=}\mathrm{-}\mathrm{23.04}\mathrm{+}\mathrm{0.48}\log {\mathit{E}}_{\mathrm{iso}}\mathrm{\pm }\mathrm{0.23}\mathit{,}$(11)where E_iso is in erg and ε_peak is in keV.

• 6)

  According to Band (2006), detectable events must have a photon flux above the threshold value C_T of the detector. Band (2006) computed this threshold for different experiments, in particular for Swift’s BAT, whose data we used here. The sensitivity was computed in the 1−1000 keV band as a function of the peak energy at the observer’s frame (${\mathit{\epsilon }}_{\mathrm{peak}}^{\mathrm{\prime }}$), taking into account the photon distribution in the considered energy interval. To facilitate the computations, the sensitivity ${\mathit{C}}_{\mathrm{T}}\mathrm{\left(}{\mathit{\epsilon }}_{\mathrm{peak}}^{\mathrm{\prime }}\mathrm{\right)}$ was fitted by the polynomial $\begin{array}{ccc}\log {\mathit{C}}_{\mathit{T}}\mathrm{\left(}{\mathit{\epsilon }}_{\mathrm{peak}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\mathrm{1.514}\mathrm{-}\mathrm{0.8801}\mathit{X}\mathrm{-}\mathrm{0.06578}{\mathit{X}}^{\mathrm{2}}& & \\ \mathrm{+}\mathrm{0.1319}{\mathit{X}}^{\mathrm{3}}\mathrm{-}\mathrm{0.02045}{\mathit{X}}^{\mathrm{4}}\mathit{,}& & \end{array}$(12)where $\mathit{X}\mathrm{=}\log {\mathit{\epsilon }}_{\mathrm{peak}}^{\mathrm{\prime }}$. In this equation the peak energy is given in keV and the photon flux C_T is given in ph cm^-2 s^-1.

To verify whether a given simulated burst triggers the detector, its average photon flux in the 1−1000 keV band was estimated by the following procedure. First, a sample including 133 LGRB with measured T₉₀ (duration in which 90% of the total burst fluence is detected) was prepared. After correcting for the rest frame, the resulting distribution $\mathit{\psi }\mathrm{\left(}{\mathit{T}}_{\mathrm{90}}^{\mathrm{0}}\mathrm{\right)}$ was fitted by a log-normal distribution characterized by a median $\log {\mathit{T}}_{\mathrm{90}}^{\mathrm{0}}\mathrm{=}\mathrm{1.213}$ and a dispersion ${\mathit{\sigma }}_{\log {\mathit{T}}_{\mathrm{90}}^{\mathrm{0}}}\mathrm{=}\mathrm{0.566}$. The distribution $\mathit{\psi }\mathrm{\left(}{\mathit{T}}_{\mathrm{90}}^{\mathrm{0}}\mathrm{\right)}$ permits one to assign the duration (in the rest frame) of the considered event, via a Monte Carlo procedure. In this case, the expected photon rate R_γ corresponding to this simulated event in the observer’s frame is ${\mathit{R}}_{\mathit{\gamma }}\mathrm{\left(}{\mathit{\epsilon }}_{\mathrm{peak}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{E}}_{\mathrm{J}}\mathit{S}\mathrm{\left(}{\mathit{\epsilon }}_{\mathrm{1}}\mathit{,}{\mathit{\epsilon }}_{\mathrm{2}}\mathrm{\right)}}{\mathrm{4}\mathit{\pi }\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathit{\theta }\mathrm{)}{\mathit{D}}_{\mathit{L}}^{\mathrm{2}}{\mathit{T}}_{\mathrm{90}}^{\mathrm{0}}}\mathrm{·}$(13)Note that in this equation there is no (1 + z) term in the denominator since it cancels with a similar term in the numerator due to the jet energy, which is also referred to the rest frame. The term S(ε₁,ε₂) is computed from the equation $\mathit{S}\mathrm{\left(}{\mathit{\epsilon }}_{\mathrm{1}}\mathit{,}{\mathit{\epsilon }}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}^{\mathrm{1000}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\mathit{B}\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}\mathrm{d}\mathit{\epsilon }\mathit{/}{\mathrm{\int }}_{\mathrm{1}}^{{\mathrm{10}}^{\mathrm{4}}}\mathit{\epsilon B}\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}\mathrm{d}\mathit{\epsilon ,}$(14)which represents the fraction of the total photons released in the event and detected in the range 1−1000 keV. Note that the limits of the integral are corrected for the rest frame and the total energy released in the event is assumed to be restricted to the normalization interval 1 keV−10 MeV. Under these conditions, detectable events satisfy the condition R_γ(ε_peak) ≥ C_T(ε_peak) and are saved as such. It is worth mentioning that although our detection condition refers to the average flux of the event, the BAT detector triggers the burst with respect to the peak flux, which is higher. In practice, this affects the estimate of the visibility function, as we see below. To compute peak fluxes a model for the burst time profile is required, which is currently not included in our code. The diversity of time profiles makes this modeling difficult, but this problem is under investigation and will be the subject of a future paper.

3. Simulations and data analysis

For each pair of the parameters characterizing the log-normal distributions of the jet energy and of the opening half-angle of the jet, that is the median and the dispersion, a series of runs were performed until the number of detected events was equal to 10⁵. The total number of runs required to satisfy this condition depends on the adopted emission model. For isotropic emission, the total number is of the order on 2 × 10⁶, while in the case of beamed emission, the number is considerably higher, on order of 7 × 10⁸.

The detected events constitute a catalog, including parameters as the redshift z, the jet energy E_J, the opening half-angle θ of the jet and the peak energy ε_peak, characterizing the hardness (or softness) of the photon spectrum. All these parameters permit one to compute the expected fluences in the Swift energy bands 15−25 keV, 25−50 keV, 50−100 keV, 100−150 keV, and 15−150 keV by the equation $\mathit{f}\mathrm{\left(}{\mathit{\epsilon }}_{\mathrm{1}}\mathit{,}{\mathit{\epsilon }}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}{\mathit{E}}_{\mathrm{J}}}{\mathrm{4}\mathit{\pi }\mathrm{(}\mathrm{1}\mathrm{-}\cos \mathit{\theta }\mathrm{)}{\mathit{D}}_{\mathit{L}}^{\mathrm{2}}}\mathit{\kappa }\mathrm{\left(}{\mathit{\epsilon }}_{\mathrm{1}}\mathit{,}{\mathit{\epsilon }}_{\mathrm{2}}\mathrm{\right)}\mathit{,}$(15)where the fraction of the bolometric energy in a given band (ε₁ − ε₂) is $\mathit{\kappa }\mathrm{\left(}{\mathit{\epsilon }}_{\mathrm{1}}\mathit{,}{\mathit{\epsilon }}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{{\mathit{\epsilon }}_{\mathrm{1}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}^{{\mathit{\epsilon }}_{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\mathit{\epsilon B}\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}\mathrm{d}\mathit{\epsilon }/{\mathrm{\int }}_{\mathrm{1}}^{{\mathrm{10}}^{\mathrm{4}}}\mathit{\epsilon B}\mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}\mathrm{d}\mathit{\epsilon }\mathit{.}$(16)After computing the fluences, they are distributed in logarithmic bins (Δlog f(ε_i,ε_i + 1) = 0.20) and their frequency distribution is computed. The simulated fluence frequency distribution in different energy bands is then compared with those derived from data of the second Swift-BAT catalog (Sakamoto et al. 2011). From the 475 events present in the catalog, we prepared a culled sample of 403 objects, discarding events with T₉₀ < 2 s and those without fluence data on all the considered energy bands. Of the excluded events, 49 have short duration, 19 have incomplete data, and 4 have an uncertain classification. The quality of the fit was measured by the χ² test, namely ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\sum _{\mathit{i}}{\mathit{\sigma }}_{\mathit{i}}^{-2}{\left({\mathit{\nu }}_{\mathit{i,}\mathrm{obs}}\mathrm{-}{\mathit{\nu }}_{\mathit{i,}\mathrm{cal}}\right)}^{\mathrm{2}}\mathit{,}$(17)where ν_i,obs and ν_i,cal are the observed and calculated fluence frequencies in the i-bin respectively.

In the next step, the code re-starts the process with a new set of parameters, which are allowed to vary within established intervals. The procedure is stopped when a minimum of the χ² test is found. In reality, the set of parameters that minimize the χ² is not exactly the same for all energy bands. Assuming a Band spectrum with the same exponents α and β for all simulated events certainly contributes to these differences and, to remedy this situation, a higher weight was given to the parameters determined from the widest energy band, that is 15−150 keV.

The calculations were performed by using the facilities of the computation center (SIGAMM) of the Observatoire de la Côte d’Azur. On average, about five CPU.h are necessary to optimize a pair of parameters for the isotropic model, while about ten CPU.h are required for the jet-model, including or excluding evolution of the mean energy.

4. Results

To check our procedure, the code was tested without the beaming correction in order to compare our results with those in the literature obtained under the same condition. In this case, steps (3) and (4) of the code (see Sect. 2.1) were skipped since there is no need to verify the alignment of the jet with the observer’s direction.

Fig. 1 Frequency distribution of fluences in the energy band 15−150 keV. Points correspond to data derived from the Swift-BAT catalog and the error bars indicate the bin width. The solid line corresponds to simulations assuming an isotropic emission model. The energy distribution defining the model is a log-normal distribution characterized by a median log Eiso = 50.96 (erg) and a dispersion σlog E = 0.92.

Figure 1 shows the best fit of the observed frequency distribution of fluences by simulated data in the energy band 15−150 keV. Simulated fluences were computed with a log-normal energy E_iso distribution, whose median is log E_iso = 50.96 ± 0.34 and whose dispersion is σ_log E = 0.92 ± 0.18. Errors indicate the uncertainties with respect to the fit of fluences in different energy bands. An early investigation by Jimenez et al. (2001), who have also assumed a log-normal distribution for the isotropic energy, led to a mean burst energy of log E_iso = 53.11 ± 0.20, about two orders of magnitude higher, but based on a sample of only eight objects. In Fig. 2 the observed frequency distribution of isotropic energies is shown in comparison with the true distribution derived from our simulations, characterized by the parameters just mentioned. The data corresponds to 170 events listed in the online Swift-BAT integrated parameters catalog (Butler et al. 2010). These data can be quite well fitted by a log-normal distribution with a median log E_iso = 52.73 and a dispersion σ_log E = 0.88. These results indicate that in reality, even in the context of the isotropic emission model the mean energy released by GRB events is about 60 times lower than the value derived directly from observation. This difference is due to the well-known Malmquist bias, which favors the detection of the more energetic events and shifts the energy distribution toward higher values, as can be seen in Fig. 2. In fact, when these selection effects are taken into account, the isotropic energy of bursts increases with the redshift according to the analysis by Wu et al. (2012), suggesting an evolution of the energy distribution function.

Fig. 2 Frequency distribution of isotropic energies. The red curve shows the true energy distribution derived from simulations. Points represent Swift-BAT data for 170 events and error bars indicate the bin width. The solid black curve is a log-normal fit of data.

Fig. 3 Visibilities derived from our simulations as a function of redshift and for three different models: isotropic emission, jet without energy evolution, and jet including evolution characterized by a parameter δ = 0.5. Note that these visibilities correspond to the sensitivity of the Swift-BAT detector.

One of the most important results derived from our simulations is the visibility function v(z), which measures the fraction of bursts generated in a given redshift interval that are effectively seen by the detector, in other words, events in the line of sight able to trigger the detector. The visibility function is essentially given by the ratio between bursts detected according to our rules in a given redshift interval and the total number of bursts generated in the same redshift interval. As we emphasized above, our detection criterion uses the average and not the peak flux. In this case, the visibility function is probably slightly underestimated. In Fig. 3 the visibility function is shown for three different sets of simulations: the isotropic model just discussed and jet models with and without energy evolution to be analyzed below. As expected, the visibility for the jet model without evolution is lower than that derived for the isotropic model, because not all events are aligned with the observer. When a moderate energy evolution is included (δ = 0.5), the resulting visibility is slightly lower than that of the jet model without evolution up to z ~ 1.5 and then increases because the events become more energetic, being able to trigger the detector. Beyond z ~ 7.5 this is true even if a comparison is made with the isotropic model (without evolution), as can be seen in Fig. 3. If the visibility function v(z) is known, the expected event rate can be computed by the equation $\frac{\mathrm{d}{\mathit{T}}_{\mathrm{grb}}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathit{K}}_{\mathrm{grb}}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{z}}_{\max }}\frac{{\mathit{R}}_{\mathrm{\ast }}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}\mathrm{)}}\frac{\mathit{v}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{d}\mathit{z}}{\sqrt{{\mathrm{\Omega }}_{\mathit{v}}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{m}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{\mathrm{3}}}}\mathit{,}$(18)where K_grb = 4πλ_grb(c/H₀)³ = 9.89 × 10¹¹λ_grb. The numerical value is obtained if the fraction of mass of the formed stars originating in GRB events (λ_grb) is given in ${\mathit{M}}_{\mathrm{\odot }}^{-1}$, the CSFR (R_∗(z)) is given in M_⊙ Mpc^-3 yr^-1, and the event rate dT_grb/dt is in yr^-1. Moreover, in the equation above, r(z) is the comoving distance in units of the Hubble radius given by $\mathit{r}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathit{z}}\frac{\mathit{dx}}{\sqrt{{\mathrm{\Omega }}_{\mathit{v}}\mathrm{+}{\mathrm{\Omega }}_{\mathrm{m}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{x}{\mathrm{)}}^{\mathrm{3}}}}\mathrm{·}$(19)For the considered energy distribution, the numerical value of the integral in Eq. (18) is 4.71 × 10^-3, implying a GRB rate of dT_grb/dt = 4.66 × 10⁹λ_grb yr^-1. Swift data indicate an all-sky frequency of 2.45 bursts/day (Gehrels et al. 2004) and, using the observed fraction of long GRBs, the observed frequency of these events is 2.0 bursts/day. Comparing this with our theoretical rate, this implies ${\mathit{\lambda }}_{\mathrm{grb}}\mathrm{=}\mathrm{1.57}\mathrm{\times }{\mathrm{10}}^{-7}{\mathit{M}}_{\mathrm{\odot }}^{-1}$. This value should be compared with the value derived by Salvaterra & Chincarini (2007)${\mathit{\lambda }}_{\mathrm{grb}}\mathrm{=}\mathrm{1.14}\mathrm{\times }{\mathrm{10}}^{-8}{\mathit{M}}_{\mathrm{\odot }}^{-1}$ based on similar assumptions, that is isotropic emission and no evolution on the burst energetics. Using our adopted local star formation rate of 0.0103 M_⊙ Mpc^-3 yr^-1, one obtains a local formation rate for LGRBs of R_grb = 1.6 Gpc^-3 yr^-1. Despite the common assumptions (isotropy and no evolution), Salvaterra & Chincarini (2007) adopted a different procedure to compute the local formation rate: their analysis was based on the burst luminosity function, while we considered the bolometric energy distribution. Moreover, our burst detection criterion is quite different, which may explain the difference of one order of magnitude between the two estimates of the local formation rate.

4.1. Jet models: local formation rates

Jet models were investigated under two main assumptions. In the first series of runs, the median of the jet energy distribution was kept constant, while in the second series of models the median was allowed to evolve with redshift according to Eq. (3). In this case, for a given value of the parameter δ, the code searches for the best fit of the observed fluence distribution by varying the median and the dispersion of the log-normal representing the true jet energy distribution, following the same procedure as described above. Figure 4 shows the best fit of the observed fluence frequency distribution by simulated data based on the jet model without evolution in the energy distribution. In this case, the model providing the best fit corresponds to a log-normal jet energy distribution with a median log E_J = 49.26 ± 0.18 and a dispersion σ_log EJ = 0.43 ± 0.21. Note that the distribution energy of the jet model is narrower than that derived for the isotropic model. As expected, the jet model implies an average energy released in the form of γ-rays about 54 times lower than the isotropic model. The reduction factor is still higher (about 3200) when the comparison is made with the mean isotropic value deduced directly from observations, alleviating, as expected, the energy budget of these explosive events.

For a jet model without evolution the numerical integral in Eq. (18) is 2.35 × 10^-5, leading to ${\mathit{\lambda }}_{\mathrm{grb}}\mathrm{=}\mathrm{3.14}\mathrm{\times }{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\odot }}^{-1}$ and a local GRB formation rate of 324 Gpc^-3 yr^-1, which is about 200 times higher than the rate derived under the assumption that the emission is isotropic.

Fig. 4 Fluence frequency distribution for the jet model without energy evolution. Points correspond to Swift-BAT data as in Fig. 1. The solid curve was derived from simulations using a log-normal for the jet energy distribution defined by a median log EJ = 49.26 (erg) and a dispersion σlog E = 0.43.

Fig. 5 Fluence frequency distribution for the jet model including energy evolution (δ = 0.5). Points correspond to Swift-BAT data as in Fig. 1. The solid curve was derived from simulations using a log-normal distribution for the jet energy distribution in which the median varies linearly with the redshift.

When evolutionary effects are included in the energy distribution, the best fit of the observed frequency distribution of fluences for δ = 0.5 was obtained with a log-normal distribution defined by the parameters log E₀ = 48.46 ± 0.26 and σ_log EJ = 0.39 ± 0.24, as shown in Fig. 5. In this particular model, the mean jet energy at z = 1 is log E_J = 49.18, while at z = 7 the mean value is log E_J = 50.48, that is about 20 times higher. In Fig. 6 the frequency distribution of jet energies corresponding to the detected events is shown. Note that the distribution is not symmetric, but has an extended high-energy wing because of the evolution of the energy distribution. From this distribution the mean expected value of the jet energy is log E_J = 49.48. It is worth mentioning that this value derived from simulations compares quite well with the mean jet energy resulting from the sample of GRBs with prominent jet breaks by Racusin et al. (2009), that is log E_J = 49.64 and with the mean jet energy derived by Schmidt (2001), that is log E_J = 49.69. The numerical integral in Eq. (18) is in this particular case equal to 2.62 × 10^-5, resulting in a formation fraction of GRBs of ${\mathit{\lambda }}_{\mathrm{grb}}\mathrm{=}\mathrm{2.82}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\odot }}^{-1}$ and in a local formation rate of 290 Gpc^-3 yr^-3. Thus, the introduction of a moderate evolution in the jet energy distribution does not significantly change the estimate of the local formation rate. It is interesting that Lu at al. (2012) also found a similar local formation rate (R_grb = 285 Gpc^-3 yr^-1), but their analysis is based on a possible anti-correlation between the opening angle and the redshift.

Fig. 6 Frequency distribution of energies for the jet model including evolution (δ = 0.5). The distribution corresponds to detected events only.

4.2. Redshift distribution

Previous investigations on the redshift distribution of GRBs led to the conclusion that without evolution either in the luminosity function or/and in the formation efficiency of GRBs, the number of events at high redshift is underestimated if the formation rate simply follows the CSFR. Since then, the number of GRBs with measured redshift has increased considerably, permitting more robust analyses. To compare the present data with our simulations, the redshift frequency distribution was calculated using 189 events given in the online Swift-BAT catalog of integrated parameters (Butler 2010).

Considering the difficulty to model the selection effects that affect the detection of the GRB host galaxy and the redshift measurement, they are not considered here. In this case, defining N_grb = dT_grb/dt, the expected redshift distribution resulting from our simulations is $\frac{\mathrm{d}{\mathit{\nu }}_{\mathrm{grb}}}{\mathrm{d}\mathit{z}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{grb}}}\frac{\mathrm{d}{\mathit{N}}_{\mathrm{grb}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathrm{d}\mathit{z}}\mathrm{·}$(20)In Fig. 7 we show the simulated redshift distributions for the three cases considered in the present investigation and data for 189 events. As expected, the isotropic and the jet models without evolution underestimate the frequency of events for z > 2. Including a moderate evolution in the jet energy distribution permits a better representation of data, in agreement with previous investigations. Note that this evolutionary jet model predicts a fraction of events for z ≥ 6 equal to 0.022, which corresponds to four events for a sample of 189 objects, in good agreement with the fact that four events with z ≥ 6 are present in the considered sample.

Fig. 7 Frequency distribution of LGRB redshifts. Data are from online Swift-BAT integrated parameter and error bars indicate the bin width. The different curves correspond to distributions derived from the present simulations: isotropic (blue curve) and jet (black curve) models without evolution and jet model including evolution (red curve).

4.3. Supernova connection

Observations seem to suggest that at least a fraction of LGRBs are linked to type Ibc supernova, whose spectra are dominated by broad absorption lines. Most of the associated bursts have isotropic energies of about 10⁴⁹ erg and are assumed to have jets with wide opening angles (E_J ~ E_iso). It is not clear yet if these low luminosity bursts (LL GRBs) constitute a particular subclass of the LGRBs (Liang et al. 2007; Guetta & Della Valle 2007). According to estimates by Soderberg et al. (2006a), the rate of LL GRBs is about 230 Gpc^-3 yr^-1, similar to the rate currently derived for LGRBs. However, Guetta & Della Valle (2007) concluded that energetic LGRBs represent only a small fraction of LL GRBs. This clearly is an open question.

The fraction by mass λ_grb of formed stars giving origin to LGRBs was estimated in the present study to be ${\mathit{\lambda }}_{\mathrm{grb}}\mathrm{=}\mathrm{2.8}\mathrm{\times }{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\odot }}^{-1}$. If the initial mass function (IMF) ζ(m) is known, λ_grb fixes the lower mass limit m_grb of progenitors of LGRBs by the relation ${\mathit{\lambda }}_{\mathrm{grb}}\mathrm{=}{\mathrm{\int }}_{{\mathit{m}}_{\mathrm{grb}}}^{{\mathit{m}}_{\mathit{u}}}\mathit{\zeta }\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{d}\mathit{m}\mathit{.}$(21)Using the IMF by Miller & Scalo (1979) normalized in the mass interval 0.1−125 M_⊙, from the equation above, one obtains that the progenitors of LGRBs must have a minimum mass of about 90 M_⊙. On the other hand, Georgy et al. (2009) estimated based on rotating stellar models computed by Meynet & Maeder (2003, 2005) that for solar metallicities the minimum mass of the progenitors of type Ic supernovae is about 39 M_⊙. Thus, from Eq. (21) one obtains ${\mathit{\lambda }}_{\mathrm{Ic}}\mathrm{\simeq }\mathrm{3.1}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{-4}{\mathit{M}}_{\mathrm{\odot }}^{-1}$ for the fraction by mass of formed stars that become type Ic supernovae. These values imply that only ~9.0% of SNIc could be associated with LGRBs, namely, only those with a very massive progenitor. Instead of using lower mass limits derived theoretically, it is possible to directly compare the present derived formation rate of LGRBs with the SNIbc rate estimated from observations. Using data by Cappellaro et al. (1999) for rates of core-collapsed supernovae and the ratio between different supernova types derived by Prieto et al. (2008), we estimated that R_SNIbc ~ 2.7 × 10⁴ Gpc^-3 yr^-1. Comparing this with the LGRB rate derived in this work, one concludes that only ~1% of type Ic supernovae could be related with LGRBs. Note that the SNIbc rate estimated from observations implies progenitor masses lower than those derived theoretically, that is around 16 M_⊙. However, these results indicate that the majority of SNIbc are not related to LGRBs, in agreement with the conclusions by Soderberg et al. (2006b). However, it is worth mentioning that recent radio observations of two SNIbc supernovae (Soderberg et al. 2010) and 2007gr (Paragi et al. 2010) not associated with a γ-flash indicate the presence of mildly relativistic outflows, which require a central engine able to power the observed flows.

We emphasize that if the progenitors of LGRBs are not single stars, our analysis is invalid. This is also true if the IMF evolves as assumed by Wang & Dai (2011), who proposed an IMF that becomes increasingly top-heavy at high reshift, producing higher LGRB rates at high z.

5. Conclusions

We reported results based on Monte Carlo simulations aimed at studying the properties of LGRBs. Mock catalogs of LGRBs derived from these simulations permitted us to predict observed fluences, and from the comparison with actual data, we were able to estimate the parameters that define the energy distribution of LGRBs. These simulations permitted us also to compute the visibility function that measures the probability for a given burst generated at a given redshift z to be able to trigger the Swift-BAT detector.

When an isotropic emission model without evolution is considered, the observed fluence distribution derived from the Swift-BAT catalog can be reproduced by an energy distribution modeled by a log-normal characterized by a median log   E_iso = 50.95 (erg) and a dispersion σ_log E = 0.92. When a comparison is made with the mean energy derived directly from observations, we realize that due to the Malmquist bias, the aforementioned intrinsic (simulated) value is about 60 times lower. The local formation rate of LGRBs derived from our simulations under the assumption of isotropic emission is R_grb = 1.6 Gpc^-3 yr^-1, an intermediate value when compared with those derived by Salvaterra & Chincarini (2007) and Pélangeon et al. (2008). The present local formation rate based on the isotropic emission model agrees with the result by Wanderman & Piran (2010), but we emphasize that in their approach evolution effects are present in the luminosity function of GRBs. The formation rate of subluminous LGRBs (under the assumption of isotropic emission) could be considerably higher. Howell & Coward (2013) derived a formation rate of ~150 Gpc^-3 yr^-1 for these events.

For beamed emission model without evolution of the energy distribution is considered, the observed fluence distribution can be fitted by simulated data when a log-normal defined by a median log E_J = 49.26 (in erg) and a dispersion σ_log   E = 0.43 represent the intrinsic energy distribution of LGRBs. As expected, these results considerably alleviate the energy requirements for GRBs models. The resulting local formation rate for this model is R_rgb = 324 Gpc^-3 yr^-1. Frail et al. (2001) obtained a local formation rate in the context of the jet model of R_grb ~ 250 Gpc^-3 yr^-1, similar to our result this agreement is fortuitous for the following reason: they took the isotropic rate given by Schmidt (2001) and corrected this value by the harmonic mean of the beaming fraction. A similar reasoning was given by Lamb et al. (2005) in their estimates of the local formation rate.

The redshift distribution predicted by the isotropic and the jet models without evolutionary effects underestimates the number of LGRBs occurring at z > 2, in agreement with previous investigations concluding the necessity of a moderate/strong evolution of the luminosity function. The model prediction of the observed redshift distribution (currently including 189 objects) can be considerably improved if the jet energy distribution evolves in the sense that at high redshift bursts are more energetic. Such an evolution was modeled in our simulations by assuming that the median of the log-normal describing the energy distribution varies as E_med,J ∝ e^δ(1 + z). Different set of simulations were performed by varying the slope dlgE_med/d(1 + z) = δ, and the best representation was found for δ = 0.5 (see Fig. 7). The ratio between the mean jet energies for events produced at z = 7 and z = 1 is about 20, indicating a moderate evolution in the energy distribution. For this model, the derived local formation is R_grb = 290 Gpc^-3 yr^-1 and the expected average jet energy is 3.0 × 10⁴⁹ erg, in agreement with the value derived by Lu et al. (2012). Nevertheless, these authors have considered the opening angle to be related to the redshift, while in our simulations both quantities were statistically independent.

The aforementioned evolutionary jet model led to a fraction by mass of formed stars originating in LGRBs of ${\mathit{\lambda }}_{\mathrm{grb}}\mathrm{=}\mathrm{2.8}\mathrm{\times }{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\odot }}^{-1}$. If the IMF does not evolve and a Miller & Scalo (1979) IMF is adopted, the derived value of λ_grb implies that the minimum mass to produce an LGRB is 90 M_⊙, indicating that only very massive stars are associated to these events. Moreover, the ratio between LGRBs and the SNIbc formation rates is in the range 0.01−0.09, supporting the idea that stars less massive than the above limit may produce an SNIbc event, but not necessarily an LGRB.

Acknowledgments

C.K. acknowledges the financial support from the University of Nice-Sophia Antipolis for her Ph.D. project. The authors also thank J.-L. Atteia and T. Regimbau for the critical reading of the manuscript.

References

All Figures

Fig. 1 Frequency distribution of fluences in the energy band 15−150 keV. Points correspond to data derived from the Swift-BAT catalog and the error bars indicate the bin width. The solid line corresponds to simulations assuming an isotropic emission model. The energy distribution defining the model is a log-normal distribution characterized by a median log Eiso = 50.96 (erg) and a dispersion σlog E = 0.92.

In the text

	Fig. 2 Frequency distribution of isotropic energies. The red curve shows the true energy distribution derived from simulations. Points represent Swift-BAT data for 170 events and error bars indicate the bin width. The solid black curve is a log-normal fit of data.
In the text

	Fig. 3 Visibilities derived from our simulations as a function of redshift and for three different models: isotropic emission, jet without energy evolution, and jet including evolution characterized by a parameter δ = 0.5. Note that these visibilities correspond to the sensitivity of the Swift-BAT detector.
In the text

	Fig. 4 Fluence frequency distribution for the jet model without energy evolution. Points correspond to Swift-BAT data as in Fig. 1. The solid curve was derived from simulations using a log-normal for the jet energy distribution defined by a median log EJ = 49.26 (erg) and a dispersion σlog E = 0.43.
In the text

	Fig. 5 Fluence frequency distribution for the jet model including energy evolution (δ = 0.5). Points correspond to Swift-BAT data as in Fig. 1. The solid curve was derived from simulations using a log-normal distribution for the jet energy distribution in which the median varies linearly with the redshift.
In the text

	Fig. 6 Frequency distribution of energies for the jet model including evolution (δ = 0.5). The distribution corresponds to detected events only.
In the text

	Fig. 7 Frequency distribution of LGRB redshifts. Data are from online Swift-BAT integrated parameter and error bars indicate the bin width. The different curves correspond to distributions derived from the present simulations: isotropic (blue curve) and jet (black curve) models without evolution and jet model including evolution (red curve).
In the text