The Luminosity Function and Formation Rate of A Complete Sample of Long Gamma-Ray Bursts

TL;DR

This study analyzes the luminosity function and formation rate of long gamma-ray bursts using a complete, redshift-known sample, revealing evolution patterns and rate densities consistent with or exceeding cosmic star formation rates.

Contribution

First application of a maximum likelihood method to a complete GRB sample with known redshifts, providing new insights into GRB evolution and rate density.

Findings

Strong redshift evolution in luminosity or density is needed to fit observations.

GRB formation rate rises steeply up to z=2 and declines thereafter.

Local GRB rate estimated at approximately 1.5 Gpc^{-3} yr^{-1}.

Abstract

We study the luminosity function and formation rate of long gamma-ray bursts (GRBs) by using a maximum likelihood method. This is the first time this method is applied to a well-defined sample of GRBs that is complete in redshift. The sample is composed of 99 bursts detected by the $Swift$ satellite, 81 of them with measured redshift and luminosity for a completeness level of $82\%$. We confirm that a strong redshift evolution in luminosity (with an evolution index of $\delta=2.22^{+0.32}_{-0.31}$) or in density ($\delta=1.92^{+0.20}_{-0.21}$) is needed in order to reproduce the observations well. But since the predicted redshift and luminosity distributions in the two scenarios are very similar, it is difficult to distinguish between these two kinds of evolutions only on the basis of the current sample. Furthermore, we also consider an empirical density case in which the GRB rate density is directly described as a broken power-law function and the luminosity function is taken to be non-evolving. In this case, we find that the GRB formation rate rises like $(1+z)^{3.85^{+0.48}_{-0.45}}$ for $z\leq2$ and is proportional to $(1+z)^{-1.07^{+0.98}_{-1.12}}$ for $z\geq2$. The local GRB rate is $1.49^{+0.63}_{-0.64}$ Gpc$^{-3}$ yr$^{-1}$. The GRB rate may be consistent with the cosmic star formation rate (SFR) at $z\leq2$, but shows an enhancement compared to the SFR at $z\geq2$.

Tables

Table 1. Table 1: Best-fitting parameters in different models.

Model	Evolution parameter	$\eta$	$a$	$b$	$\log L_c$	${\chi }_{\mathrm{total}}^2$	AIC
		(${10}^{-8}$ $M_{\odot }^{-1}$)			(erg $s^{-1}$)
No evolution	$\cdots$	${9.00}_{-2.67}^{+2.22}$	$-{0.36}_{-0.16}^{+0.16}$	$-{1.28}_{-0.09}^{+0.09}$	${52.17}_{-0.14}^{+0.18}$	20169.7	20177.7
Luminosity evolution	$\delta ={2.22}_{-0.31}^{+0.32}$	${7.35}_{-2.06}^{+2.11}$	$-{0.72}_{-0.05}^{+0.05}$	$-{1.42}_{-0.24}^{+0.24}$	${52.12}_{-0.26}^{+0.26}$	20086.1	20096.1
Density evolution	$\delta ={1.92}_{-0.21}^{+0.20}$	${4.60}_{-1.63}^{+1.67}$	$-{0.69}_{-0.06}^{+0.05}$	$-{1.75}_{-0.29}^{+0.31}$	${53.35}_{-0.15}^{+0.16}$	20083.2	20093.2
Empirical density	${\psi }_{\mathrm{GRB}}\left(0\right)={1.49}_{-0.64}^{+0.63}$ $n_1={3.85}_{-0.45}^{+0.48}$ $n_2=-{1.07}_{-1.12}^{+0.98}$ $z_c={2.33}_{-0.24}^{+0.39}$	$\cdots$	$-{0.69}_{-0.06}^{+0.06}$	$-{1.76}_{-0.33}^{+0.34}$	${53.32}_{-0.17}^{+0.18}$	20079.0	20093.0

Equations

$$\mathcal{L}=\exp (-N_{\exp })\prod _{i=1}^{N_{\mathrm{obs}}}\Phi (L_i,z_i,t_i),$$

$$\Phi \left(L,z,t\right)=\frac{d^{3}N}{dtdzdL}=\frac{\Delta \Omega }{4\pi }\frac{\psi (z)}{1+z}\frac{dV(z)}{dz}\varphi (L),$$

$${\chi }_{\mathrm{Swift}}^2=-2\sum _{i=1}^{N_{\mathrm{obs}}}\ln \left[\Phi (L_i,z_i,t_i)\right]+2N_{\exp }.$$

$$N_{\exp }=\frac{\Delta \Omega T}{4\pi }{\int }_0^{z_{\max }}\frac{\psi (z)}{1+z}\frac{dV(z)}{dz}dz{\int }_{\max [L_{\min },L_{\lim }(z)]}^{L_{\max }}\varphi (L)dL,$$

$$L_{\lim }(z)=4\pi D_L^2(z)P_{\lim }\frac{{\int }_{1/(1+z)\mathrm{keV}}^{10^4/(1+z)\mathrm{keV}}EN(E)dE}{{\int }_{15\mathrm{keV}}^{150\mathrm{keV}}N(E)dE},$$

$$N^{\exp }(P_1<P<P_2)=0.5\times 10.8$$

$$\times {\int }_{L(P_1,z)}^{L(P_2,z)}\varphi (L)dL.$$

$${\chi }_{\mathrm{Fermi}}^2=\sum _i^n\frac{{\left(N_i^{\mathrm{obs}}-N_i^{\exp }\right)}^2}{{\sigma }_i^2},$$

$${\chi }_{\mathrm{total}}^2={\chi }_{\mathrm{Swift}}^2+{\chi }_{\mathrm{Fermi}}^2.$$

$$\phi(L)=\frac{A}{\ln(10)L}\left\{\begin{array}[]{l}{\left(\frac{{{L}}}{{{L_{c}}}}\right)^{a}};\,\,{L}\leq{L_{c}}\\ {\left(\frac{{{L}}}{{{L_{c}}}}\right)^{b}};\,\,{L}>{L_{c}}\;,\end{array}\right.$$

$$P({\mathcal{M}}_{\alpha })=\frac{\exp (-\mathrm{AI}{\mathrm{C}}_{\alpha }/2)}{\exp (-\mathrm{AI}{\mathrm{C}}_1/2)+\cdot \cdot \cdot +\exp (-\mathrm{AI}{\mathrm{C}}_{\mathrm{N}}/2)}.$$

$${\psi }_{\star }(z)=\frac{0.0157+0.118z}{1+(z/3.23{)}^{4.66}}.$$

$$\psi(z)={\psi_{\rm GRB}}(0)\left\{\begin{array}[]{l}{(1+z)^{{n_{1}}}};\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;z\leq{z_{c}}\\ {(1+{z_{c}})^{{n_{1}}-{n_{2}}}}{(1+z)^{{n_{2}}}};\;z>{z_{c}}\;,\end{array}\right.$$

Full text

††pubyear: 2018††pagerange: The Luminosity Function and Formation Rate of A Complete Sample of Long Gamma-Ray Bursts–The Luminosity Function and Formation Rate of A Complete Sample of Long Gamma-Ray Bursts

1 Introduction

Gamma-ray bursts (GRBs) are the most energetic explosions in the universe, which can be detected up to extremely high redshifts. In theory, long GRBs with durations $T_{\rm 90}>2$ s (where $T_{\rm 90}$ is the time interval observed to contain $90\%$ of the prompt emission; Kouveliotou et al. 1993) are believed to originate from the core collapse of massive stars (e.g., Woosley 1993; Paczyński 1998; Woosley & Bloom 2006), an idea given significant support from some confirmed associations between long GRBs and supernovae (e.g., Hjorth et al. 2003; Stanek et al. 2003). This collapsar model implies that the GRB formation rate should in principle trace the cosmic star formation rate (SFR; Totani 1997; Wijers et al. 1998; Lamb & Reichart 2000; Porciani & Madau 2001; Piran 2004; Zhang & Mészáros 2004; Zhang 2007). However, the $Swift$ observations seem to indicate that the GRB rate does not closely follow the SFR but is actually enhanced by some unknown mechanisms at high-$z$ (Daigne et al., 2006; Guetta & Piran, 2007; Le & Dermer, 2007; Salvaterra & Chincarini, 2007; Kistler et al., 2008, 2009; Li, 2008; Yüksel et al., 2008; Salvaterra et al., 2009, 2012; Campisi et al., 2010; Qin et al., 2010; Wanderman & Piran, 2010; Cao et al., 2011; Virgili et al., 2011; Elliott et al., 2012; Lu et al., 2012; Robertson & Ellis, 2012; Tan et al., 2013; Wang, 2013; Wei et al., 2014; Tan & Wang, 2015; Deng et al., 2016; Wei & Wu, 2017; Paul, 2018).111Using the $C^{-}$ statistical method proposed by Lynden-Bell (1971), Pescalli et al. (2015) and Yu et al. (2015) found a relative excess of the GRB formation rate with respect to the SFR at low redshifts. But then Pescalli et al. (2016) showed that if the $C^{-}$ method is applied to incomplete GRB samples it can misleadingly lead to an excess of the GRB rate at $z\leq 1$. In addition, some works performed spectro-photometric studies on the properties (stellar mass, SFR, and metallicity) of long GRB host galaxies of different complete GRB samples and compared them to the ones of typical star-forming galaxies selected by galaxy surveys (e.g., Vergani et al. 2015; Japelj et al. 2016; Perley et al. 2016; Palmerio et al. 2019). All their results clearly suggested that at $z\leq 1$ only a small fraction of the star formation produces GRBs. Several evolution models have been proposed to explain the observed enhancement, such as the GRB rate density evolution (Kistler et al., 2008, 2009), cosmic metallicity evolution (Langer & Norman, 2006; Li, 2008), and an evolution in the GRB luminosity function (Virgili et al., 2011; Salvaterra et al., 2012; Tan et al., 2013; Tan & Wang, 2015; Paul, 2018).

However, it should be emphasized that our knowledge about the properties of long GRBs and their evolution with cosmic time is still hindered by the fact that most of the observed $Swift$ GRBs are without redshift. Indeed, only $\sim 1/3$ of all $Swift$ GRBs have redshift determinations. Most of previous researches adopted incomplete redshift samples to derive the GRB luminosity function and redshift distribution. Given the low completeness level in redshift measurement, the possible observational biases may have remarkable effect on sharping the GRB redshift distributions (Fiore et al., 2007; Salvaterra et al., 2012). Therefore, it is necessary to consider using an unbiased complete sample of long GRBs that is capable of adequately representing this class of object to study their distributions through cosmic time. Salvaterra et al. (2012) defined a complete flux-limited sample of $Swift$ long GRBs which, despite containing a relatively small sample size, has a completeness level in redshift determination of 90%. The high level of redshift completeness enabled them for the first time to constrain the GRB luminosity function and its evolution in an unbiased way. They found that either a luminosity evolution with $\delta=2.1\pm 0.6$ or a density evolution with $\delta=1.7\pm 0.5$ can well reproduce the observations. However, they can not discriminate between these two scenarios. Recently, Pescalli et al. (2016) revised the complete sample of Salvaterra et al. (2012) and then extended it with new bursts that have favorable observing conditions for redshift determination and that are bright in the 15–150 keV $Swift$/BAT band. The updated sample is composed of 99 bursts, 81 of them with known redshift and luminosity for a completeness level of 82%. Pescalli et al. (2016) adopted the Lynden-Bell $C^{-}$ method to derive the luminosity function and formation rate of GRBs from this updated complete $Swift$ sample. A strong evolution in luminosity $L(z)\propto(1+z)^{2.5}$ was found.

There are several algorithms to derive the luminosity function for a specific kind of astronomical sources. A classical approach to determine the luminosity function is based on the $1/V_{\rm MAX}$ method of Schmidt (1968) applied to redshift bins. However, it is known that this method would introduce bias if there is strong evolution within the bins. Moreover, given the relatively small number of GRBs and their wide redshift and luminosity range spanned, binning would lead to a loss of information. Other non-parametric methods (e.g., the $C^{-}$ method; Lynden-Bell 1971) usually require certain uniformity of data coverage to be applicable. In practice, the GRB data are truncated by the flux sensitivity limit of the detector. It is very difficult to parametrize the sensitivity of the detector and to construct a uniformly distributed GRB sample. The maximum-likelihood algorithms (Marshall et al., 1983) are therefore preferable for the GRB problems, because these methods are more flexible in modeling the systematical uncertainties and are less limited by the conditions of a given sample. Wanderman & Piran (2010) adopted a maximum likelihood estimator to obtain the GRB luminosity function and formation rate. They analyzed a sample of $\sim 100$ $Swift$ GRBs with known redshifts, which is possibly suffering from incompleteness.

In this work, we make use of the high completeness of the updated sample presented in Pescalli et al. (2016) to constrain the GRB luminosity function and redshift distribution. Compared with previous works using incomplete samples, which relied on the assumption that bursts lacking redshift information strictly follow the redshift distribution of bursts with measured redshifts, our present work has the advantage of being independent of this assumption. Additionally, to include systematics and unknowns in the statistical inference, we apply the maximum likelihood method, for the first time, to analyze the complete GRB sample. The rest of the paper is organized as follows. In Section 2, we describe the complete sample at our disposal. In Section 3, we illustrate the maximum likelihood method used for our analysis. Our models and analysis results are presented in Section 4. Lastly, we draw a brief summary in Section 5. Throughout this paper we adopt a standard $\Lambda$CDM cosmological model with $\Omega_{m}=0.3$, $\Omega_{\Lambda}=0.7$, and $H_{0}=70$ km s*-1* Mpc*-1*.

2 The Sample

Since the launch of the $Swift$ satellite (Gehrels et al., 2004), the number of measured GRB redshifts has increased rapidly. But we still have to face the problem that the bursts with redshift measurements account for only $\sim 1/3$ of all GRBs. The low completeness level in redshift determination will undoubtedly lead to biases in the statistical analysis of GRBs. Jakobsson et al. (2006) thus proposed some criteria to select long GRBs which have favorable observing conditions for redshift measurement. Salvaterra et al. (2012) built a complete sample of $Swift$ long GRBs (called BAT6), which is composed of 58 GRBs matching the criteria of Jakobsson et al. (2006) and having 1-s peak photon flux $P\geq 2.6$ ph cm*-2* s*-1* (integrated in the 15–150 keV BAT energy band). 52 of them have measured redshift so that the completeness level is $\sim 90$%. Pescalli et al. (2016) revised the BAT6 sample and extended it with additional bursts that satisfy its selection criteria. The BAT6 extended (BAT6ext) sample contains 99 GRBs up to 2014 July, of which 81 bursts have measured $z$ and $L$. Its completeness in redshift is $\sim 82$%. In the following, we will use the BAT6ext sample to investigative the GRB luminosity function and redshift distribution.

3 Maximum Likelihood Analysis

In order to constrain the model parameters, a maximum likelihood method first introduced by Marshall et al. (1983) is adopted. The likelihood function $\mathcal{L}$ is defined by the expression (Chiang & Mukherjee, 1998; Narumoto & Totani, 2006; Ajello et al., 2009, 2012; Abdo et al., 2010; Zeng et al., 2014, 2016)

[TABLE]

where $N_{\rm exp}$ is the expected number of GRB detections, $N_{\rm obs}$ is the number of the observed sample, and $\Phi(L,z,t)$ is the observed rate of GRBs per unit time at redshift $\in(z,z+dz)$ with luminosity $\in(L,L+dL)$, which can be expressed as

[TABLE]

where $\Delta\Omega=1.33$ sr is the solid angle covered on the sky by $Swift$, $\psi(z)$ is the comoving formation rate of GRBs in units of $\rm Gpc^{-3}$ $\rm yr^{-1}$, $(1+z)^{-1}$ is the cosmic expansion factor, and $\phi(L)$ is the normalized GRB luminosity function. And $dV(z)/dz=4\pi cD_{L}^{2}(z)/[H(z)(1+z)^{2}]$ is the comoving volume element, where $D_{L}(z)$ is the luminosity distance and $H(z)=H_{0}[\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}]^{1/2}$ is the Hubble parameter. Transforming the likelihood function to the standard expression $\chi^{2}=-2\ln L$, we obtain the $\chi^{2}$ distribution for the complete $Swift$ sample:

[TABLE]

Considering the flux threshold used to defined the BAT6ext (i.e., $P_{\rm lim}=2.6$ ph cm*-2* s*-1* in the 15–150 keV energy band), the expected number of GRBs can be estimated by

[TABLE]

where $T\sim 9.3$ yr is the observational period of $Swift$ that covers the BAT6ext sample. Since $z<10$ for the current GRB sample, we adopt a maximum GRB redshift $z_{\rm max}=10$. The luminosity function is assumed to extend between minimum and maximum luminosities $L_{\rm min}=10^{49}$ erg $\rm s^{-1}$ and $L_{\rm max}=10^{55}$ erg $\rm s^{-1}$ (Pescalli et al., 2015). The luminosity threshold appearing in Equation (4) can be calculated by

[TABLE]

where $N(E)$ is the observed photon spectrum. To describe the typical GRB spectrum, we use a Band function with low- and high-energy spectral indices $-1$ and $-2.3$, respectively (Band et al., 1993; Preece et al., 2000; Kaneko et al., 2006). The spectral peak energy $E_{p}$ is obtained through the $E_{p}$–$L$ correlation (Yonetoku et al., 2004; Nava et al., 2012): $\log\left[E_{p}(1+z)\right]=-25.33+0.53\log L$, where $L$ represents the isotropic peak luminosity.

Similar to what Salvaterra et al. (2012) did in their treatment, we optimize the model free parameters by jointly fitting the observed redshift and luminosity distributions of bursts in the BAT6ext sample and the observed differential peak-flux number counts in the 50–300 keV band of $Fermi$/GBM (Gruber et al., 2014; von Kienlin et al., 2014; Narayana Bhat et al., 2016).222 The $Fermi$/GBM GRB data are available at https://fermi.gsfc.nasa.gov/ssc/data/access/gbm/. The $Fermi$/GBM sample contains $\sim 10.8$ yr of observation (up to 2019 May) with an average exposure factor of $\sim 0.5$. To avoid the complication that would arise from the use of a detailed treatment of the $Fermi$/GBM threshold, we select 1375 $Fermi$/GBM long GRBs with peak flux $P\geq 1.0$ ph cm*-2* s*-1* in the 50–300 keV energy band. The expected number of events in each peak-flux bin $P\in(P_{1},P_{2})$ should be

[TABLE]

The $\chi^{2}$ value for the $Fermi$ sample is then given by

[TABLE]

where $n$ is the number of $P$ bins, and $N^{\rm obs}$ and $\tilde{N}^{\rm exp}$ are the observed and expected numbers of GRBs in bin $i$, respectively. For the observed number $N^{\rm obs}_{i}$ in bin $i$, the statistical error of $N^{\rm obs}_{i}$ is usually considered to be the Poisson error, i.e., $\sigma_{i}=\sqrt{N^{\rm obs}_{i}}$, which denotes the 68% Poisson confidence intervals for the binned events. Here the differential peak-flux distribution is treated as a sum of independent measurements in the different 20 $P$ bins with the same width $\Delta\log P$ in $\log(P)$ space (i.e., $\Delta\log P=[\log(P_{\rm max})-\log(P_{\rm min})]/20$). Note that taking different values for $\Delta\log P$ has little impact on the best-fitting results. It should be emphasized that the measure of the local GRB rate density is only determined by the observed total number of $Fermi$ GRBs detected above $P=1.0$ ph cm*-2* s*-1*.

Therefore, we can define a new $\chi^{2}_{\rm total}$ function that combining the maximum likelihood analysis and the constraints from the number counts of $Fermi$:

[TABLE]

The complete $Swift$ sample can provide a powerful test for the existence and the level of redshift evolution of long GRBs, while the fit to the $Fermi$ number counts enables us to infer the local GRB rate density and to better constrain the luminosity function free parameters (Salvaterra et al., 2012). For each model, we optimize the free parameters using the Markov Chain Monte Carlo (MCMC) technology,333We adopt the MCMC code from CosRayMC (Liu et al., 2012), which itself was adopted from the COSMOMC package (Lewis & Bridle, 2002). which is widely employed to give multidimensional parameter constraints from the observational data. In practice, this means we will find the parameter values that minimize $\chi^{2}_{\rm total}$, which yields the best-fitting parameters and their corresponding $1\sigma$ uncertainties. It is worth stressing that the derived best-fitting parameters also give a good fit of the $Swift$ peak-flux number counts in the 15–150 keV band.

4 Models and Analysis Results

Here we explore the expression of the GRB luminosity function in a broken power law, which is widely adopted in the literature:

[TABLE]

where $A$ is a normalization constant, $L_{c}$ is the break luminosity, and $a$ and $b$ are the faint- and bright-end power-law indices, respectively.

Using the maximum likelihood analysis method described in Section 3, we can optimize the values of each model’s free parameters, including the GRB luminosity function, the GRB formation efficiency $\eta$, and the evolution parameter. Table 1 lists the best-fitting parameters together with their $1\sigma$ confidence level for different models. In the last two columns, we report the total $\chi_{\rm total}^{2}$ value and the Akaike information criterion (AIC) score, respectively. For each fitted model, the AIC is given by ${\rm AIC}\equiv\chi_{\rm total}^{2}+2k$, where $k$ is the number of free parameters (Akaike, 1974; Liddle, 2007). If there are two or more models for the data, $\mathcal{M}_{\rm 1},\mathcal{M}_{\rm 2},...,\mathcal{M}_{\rm N}$, and they have been separately fitted, the one with the least AIC score is the one most favoured by this criterion. A more quantitative ranking of models can be calculated as follows. With ${\rm AIC}_{\alpha}$ characterizing model $\mathcal{M}_{\alpha}$, the un-normalized confidence that this model is true is the Akaike weight $\exp(-{\rm AIC_{\alpha}}/2)$. The relative probability that $\mathcal{M}_{\alpha}$ is statistically preferred is

[TABLE]

One-dimensional (1-D) probability distributions and two-dimensional (2-D) regions with the 1–2$\sigma$ contours corresponding to the parameters in different models are shown in Figures 1-4. One can see from these plots that all of the model parameters are well constrained.

4.1 No evolution model

In the first simple (no-evolution) model, we assume that the GRB formation rate purely follows the cosmic SFR, $\psi_{\star}(z)$, i.e., $\psi(z)=\eta\psi_{\star}(z)$, and that their luminosity function does not evolve with redshift, i.e., $L_{c,z}=L_{c,0}={\rm const}$. The factor $\eta$ denotes the GRB formation efficiency in units of ${\rm M}_{\odot}^{-1}$. The observed SFR $\psi_{\star}(z)$ (in units of ${\rm M}_{\odot}$ $\rm yr^{-1}$ $\rm Mpc^{-3}$) is commonly parameterized with the form (Hopkins & Beacom, 2006; Li, 2008):

[TABLE]

The 1-D probability distributions and 2-D regions with the 1–2$\sigma$ contours corresponding to four parameters in the no-evolution model are presented in Figure 1.

Figure 5 shows the $z$ and $L$ distributions of 81 GRBs in the BAT6ext complete sample. The expectation from the no-evolution case (yellow dot-dashed lines) do not provide a good representation of the observed $z$ and $L$ distributions of the BAT6ext sample. Particularly, the rate of GRBs at high redshift is clearly under-predicted and the reproduce of the $L$ distribution is not as good as those of the luminosity evolution model or density evolution model, more fully described below. According to the AIC model selection criterion, we can safely discard this model as having a probability of only $\sim 10^{-19}$ of being correct compared to the other three models.

4.2 Luminosity evolution model

This model assumes that an evolution in the GRB luminosity function can enhance the number of GRB detections at high-z. In this case, while the GRB formation rate is still proportional to the cosmic SFR, the break luminosity in the GRB luminosity function increases with redshift as $L_{c,z}=L_{c,0}(1+z)^{\delta}$. In Figure 2, we also display the 1-D probability distributions and 1–2$\sigma$ constraint contours for five parameters in this model. We find that a strong luminosity evolution with $\delta=2.22^{+0.32}_{-0.31}$ is required to reproduce both the observed $z$ and $L$ distributions of 81 bursts in the BAT6ext sample (blue dashed lines in Figure 5). Based on the same sample, Pescalli et al. (2016) found a strong evolution in luminosity ($\delta=2.5$) through the Lynden-Bell $C^{-}$ method, in good agreement with our results from the maximum likelihood analysis method. Using the AIC model selection criterion, we find that this model is somewhat disfavored statistically compared to the other three models, with a relative probability of $\sim 10\%$.

4.3 Density evolution model

This model assumes that an evolution of the GRB formation rate can also provide an enhancement of the high-z GRB detection. In this scenario, while the break luminosity in the GRB luminosity function is still a constant, the GRB rate follows the cosmic SFR in conjunction with an additional evolution characterized by $(1+z)^{\delta}$, i.e., $\psi(z)=\eta\psi_{\star}(z)(1+z)^{\delta}$. The best-fitting parameters and their constraint contours are shown in Figure 3. We find that a strong density evolution with $\delta=1.92^{+0.20}_{-0.21}$ reproduces the observed $z$ and $L$ distributions (red solid lines in Figure 5) quite well. On the basis of the AIC model selection criterion, we find that among four different models, this one is statistically preferred with a relative probability $\sim 43\%$.

The large value of $\delta$ means an obvious shift of the peak of the GRB formation rate toward a higher redshift with respect to stars. In this following section, we will further investigate this issue by directly adopting an empirical function as the GRB rate density, without making any assumption on the relation between the GRB rate and the SFR. Later on we will compare this empirical function with models for the GRB rate that follow the SFR.

4.4 Empirical density model

We approximate the expression of the GRB formation rate with an empirical broken power-law function (Wanderman & Piran, 2010):

[TABLE]

where $\psi_{\rm GRB}(0)$ is the local GRB formation rate and $z_{c}$ is the break redshift. Same as the density evolution model, there is no evolution of the GRB luminosity function in this case. Figure 4 displays the constraint results on the parameters $\psi_{\rm GRB}(0)$, $a$, $b$, $L_{c,0}$, $z_{c}$, $n_{1}$, and $n_{2}$. We find that the expectation from the empirical density case also provide a good representation of the observed $z$ and $L$ distributions (green dotted lines in Figure 5). According to the AIC, the empirical density model is slightly favored compared to the density evolution model, but the differences are statistically insignificant ($\sim 47\%$ for the former versus $\sim 43\%$ for the latter). In both of these two models, the GRB formation rate is found to peak at a higher redshift with respect to the no-evolution model. The intrinsic GRB formation rate as a function of redshift in different models are displayed in Figure 6. The comparison between the no-evolution and empirical density models shows that the GRB rate may follow the SFR for $z\la 2$, but shows an enhancement compared to the SFR for $z\ga 2$.

5 Conclusion and Discussion

In this work, we try to investigate the properties of long GRBs and their evolution with redshift. To achieve this aim, a maximum likelihood method (Marshall et al., 1983) is adopted, whose specific version has already been applied to GRBs (Wanderman & Piran, 2010). Here we apply this method, for the first time, to a carefully selected sample of long GRBs that is complete in peak flux and 82% complete in redshift. This sample is composed of GRBs detected by the $Swift$ satellite with favorable observing conditions for redshift determination and with peak photon fluxes $P\geq 2.6$ ph cm*-2* s*-1* (Pescalli et al., 2016). It contains 99 bursts with a completeness of $\sim 82\%$ (81 out of 99 bursts with known $z$ and $L$).

Based on the complete GRB sample, we directly construct the GRB luminosity function and redshift distribution in the frameworks of different evolution models using the maximum likelihood method that performs an analysis for the redshift and luminosity distributions of the GRB sample and the MCMC technology that provides the best-fitting parameters (see Table 1) and the probability density distributions of the parameters in each model (see Figures 1-4). According to the AIC model selection criterion, we confirm that the no-evolution model can be safely excluded. That is, GRBs must have experienced some kind of evolution to become more luminous or more population in the past than present day. In order to account for the observed distributions, the GRB luminosity should increase to $(1+z)^{2.22^{+0.32}_{-0.31}}$ or the GRB rate density should increase to $(1+z)^{1.92^{+0.20}_{-0.21}}$ with respect to the known cosmic SFR. But since the luminosity and density evolution scenarios predict very similar distributions, we can not distinguish between these two kinds of evolutions simply on the basis of the sample used in this study. These results are in good agreement with those of other works (Salvaterra et al., 2012; Pescalli et al., 2016).

Note that both the luminosity and density evolution models explored here are based on the assumption that the GRB formation rate is related to a given SFR. With different SFR models, the constraint results may change in some degree (see Virgili et al. 2011; Hao & Yuan 2013; Wei et al. 2016). Therefore, we also consider an empirical density case in which the GRB rate density is approximated as a broken power-law function, rather than being related to the SFR. In this case, we find that the GRB rate rapidly increases at $z\la 2$ ($n_{1}\simeq 3.85$) and then shows slowly decreases at $z\ga 2$ ($n_{2}\simeq-1.07$). The GRB rate is compatible, of course, with a constant rate at $z\ga 2$. The local formation rate of GRBs is $1.49^{+0.63}_{-0.64}$ Gpc*-3* yr*-1*, which is consistent with the result of Wanderman & Piran (2010). By comparing the derived intrinsic redshift distributions in the no-evolution and empirical density models, we find that while the GRB rate may be consistent with the SFR at $z\la 2$, its high-redshift slope is shallower than the steep decline of the SFR at $z\ga 2$.

The total (complete) sample presented by Pescalli et al. (2016) comprises 99 GRBs. However, only 81 of them with known redshift and luminosity have been used in our analysis. To investigate the impact of the inclusion of the other 18 bursts without known redshift on our conclusions, we apply the luminosity correlation (Nava et al., 2012), $\log\left[E_{p}(1+z)\right]=-25.33+0.53\log L$, to estimate their pseudo redshifts. We use the observed peak flux and $E_{p}$ of the 18 bursts to calculate the rest-frame peak energies and the isotropic peak luminosities for different redshifts. By requiring the bursts enter the $2\sigma$ (or $1\sigma$) region of the correlation, we derive the lower limits of redshifts and then the corresponding luminosities for these 18 bursts. In Figure 7, we present the $z$ and $L$ distributions of 81 GRBs with measured redshifts (red dots) and of 99 GRBs (including 81 bursts with measured redshifts and 18 ones with pseudo redshifts), respectively. Blue squares and green diamonds correspond to the cases of requiring the 18 GRBs enter the $1\sigma$ and $2\sigma$ regions of the luminosity correlation, respectively. One can see from this plot that the $z$ and $L$ distributions of these three cases are almost the same, we can therefore conclude that the inclusion of this 18% of bursts without known $z$ would not change the main conclusions of our paper.

Acknowledgements

We are grateful to the anonymous referee for insightful comments. We also thank Bin-Bin Zhang for helpful discussion on the average exposure factor of $Fermi$/GBM. This work is partially supported by the National Natural Science Foundation of China (grant Nos. 11603076, 11673068, 11725314, 11703094, and U1831122), the Youth Innovation Promotion Association (2017366), the Key Research Program of Frontier Sciences (QYZDB-SSW-SYS005), the Strategic Priority Research Program “Multi-waveband gravitational wave Universe” (grant No. XDB23000000) of Chinese Academy of Sciences, and the “333 Project” and the Natural Science Foundation (grant No. BK20161096) of Jiangsu Province.