EXTRAGALACTIC HIGH-ENERGY TRANSIENTS: EVENT RATE DENSITIES AND LUMINOSITY FUNCTIONS

Quantifying the redshift-dependent event rate density and LF of a certain type of transient is a challenging task. The observed events and their redshift and luminosity distributions are the results of the convolution of both (likely redshift-dependent) LF and intrinsic redshift distribution of the events, with the proper correction of the instrumental sensitivity threshold, field of view, and operational time. All these complications may be disentangled with a large enough observational sample. However, for most transients discussed in this paper (except GRBs), the number of observational samples is too small to perform such a task. In order to cross-compare various types of transients, in this paper we first ignore the possible redshift evolution of the LFs of all the transients and use the data to construct LFs of each type of transient. This allows us to separate LF and redshift distribution f(z). In principle, the LF could be redshift dependent. For example, for LGRBs for which we have collected a large enough sample, evidence of such an evolution effect has been collected. As a result, the LF we construct in this paper only carries the meaning in the "global" sense. We hereby define all the LFs constructed without considering redshift evolution as GLFs. In order to use the data to construct the GLF, we also assume that the events with the same luminosity share the same other properties (e.g., spectral properties and detector parameters). This makes the observed events good indicators of the underlying general population. Also, since there is a wide range of the spectral peak energy (E_peak) distribution for different transients studied in this paper, we try as much as possible to apply the k-corrected bolometric luminosity (1–10⁴ keV) using the measured spectral properties of the transients (Equation (28)). The only exceptions are the TDEs detected by ROSAT, Chandra, and XMM-Newton, whose narrow bandpass does not allow a precise inference of the global spectral parameters. For these events, we use a uniform k-correction parameter 1.4 (Equation (30)).

For a certain type of transient, we define the local specific event rate density (local event rate density per unit luminosity) as

$$\begin{eqnarray}&&{\rho }_{0,L}=\displaystyle \frac{d{\rho }_{0}}{{dL}},\end{eqnarray} \tag{ 1 }$$

where L denotes bolometric luminosity at the peak time hereafter. For a detector with flux sensitivity F_th, field of view Ω, and operational time T, the detected number of events in the luminosity interval from L to $(L+{dL})$ is

$$\begin{eqnarray}{dN} & = & \displaystyle \frac{{\rm{\Omega }}T}{4\pi }{\rho }_{0,L}{dL}{\displaystyle \int }_{0}^{{z}_{{\rm{max}}}(L)}\displaystyle \frac{f(z)}{1+z}\displaystyle \frac{{dV}(z)}{{dz}}{dz}\\ & = & \displaystyle \frac{{\rm{\Omega }}T}{4\pi }{\rho }_{0,L}g(L){dL},\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&g(L)={\displaystyle \int }_{0}^{{z}_{{\rm{max}}}(L)}\displaystyle \frac{f(z)}{1+z}\displaystyle \frac{{dV}(z)}{{dz}}{dz}.\end{eqnarray} \tag{ 3 }$$

The dimensionless function f(z) describes the redshift-dependent event rate density, i.e.,

$$\begin{eqnarray}&&{\rho }_{L}(z)={\rho }_{0,L}f(z).\end{eqnarray} \tag{ 4 }$$

The redshift-dependent specific comoving volume reads (for the standard ΛCDM cosmology)

$$\begin{eqnarray}&&\displaystyle \frac{{dV}(z)}{{dz}}=\displaystyle \frac{c}{{H}_{0}}\displaystyle \frac{4\pi {D}_{{\rm{L}}}^{2}}{{(1+z)}^{2}{[{{\rm{\Omega }}}_{M}{(1+z)}^{3}+{{\rm{\Omega }}}_{{\rm{\Lambda }}}]}^{1/2}}.\end{eqnarray} \tag{ 5 }$$

For a particular L, the maximum redshift z_max(L), which defines the maximum volume inside which an event with luminosity L can be detected, can be defined by the sensitivity threshold F_th via

$$\begin{eqnarray}&&{F}_{{\rm{th}}}=\displaystyle \frac{L}{4\pi {D}_{{\rm{L}}}^{2}({z}_{{\rm{max}}})k},\end{eqnarray} \tag{ 6 }$$

where k is a correction factor, which corrects the observed flux in the detector's energy band (e₁, e₂) to a wide band in the rest frame (e.g., 1–10⁴ keV for GRBs; see Equation (29)).

Technically, it is easier to evaluate numbers in the logarithmic luminosity bins. Equation (2) can be also written as

$$\begin{eqnarray}&&{dN}=\displaystyle \frac{{dN}}{d\mathrm{log}L}d\mathrm{log}L=\displaystyle \frac{{\rm{\Omega }}T}{4\pi }(\mathrm{ln}10){\rho }_{0,L}g(L){Ld}(\mathrm{log}L).\end{eqnarray} \tag{ 7 }$$

Suppose that ΔN events are detected in a finite logarithmic luminosity bin from $\mathrm{log}L$ to $\mathrm{log}L+{\rm{\Delta }}(\mathrm{log}L);$ one then has

$$\begin{eqnarray}&&{\rho }_{0,L}\simeq \displaystyle \frac{4\pi }{{\rm{\Omega }}T}\displaystyle \frac{1}{\mathrm{ln}10}\displaystyle \frac{1}{g(L)}\displaystyle \frac{1}{L}\displaystyle \frac{{\rm{\Delta }}N}{{\rm{\Delta }}(\mathrm{log}L)}.\end{eqnarray} \tag{ 8 }$$

The LF of a certain type of transient can be defined as

$$\begin{eqnarray}&&N(L){dL}\propto {\rm{\Phi }}(L){dL},\end{eqnarray} \tag{ 9 }$$

with the integration of Φ(L) normalized to unity, i.e.,

$$\begin{eqnarray}&&{\displaystyle \int }_{{L}_{{\rm{min}}}}^{{L}_{{\rm{max}}}}{\rm{\Phi }}(L){dL}=1,\end{eqnarray} \tag{ 10 }$$

where L_min and L_max are the minimum and maximum values of the luminosity distribution. One can define the local event rate density above a certain luminosity L, i.e.,

$$\begin{eqnarray}{\rho }_{0,\gt L} & = & {\displaystyle \int }_{L}^{{L}_{{\rm{max}}}}{\rho }_{0,L}{dL}\\ & \simeq & \displaystyle \sum _{\mathrm{log}L}^{\mathrm{log}{L}_{{\rm{max}}}}\displaystyle \frac{4\pi }{{\rm{\Omega }}T}\displaystyle \frac{1}{\mathrm{ln}10}\displaystyle \frac{1}{g(L)}\displaystyle \frac{{\rm{\Delta }}N}{{\rm{\Delta }}(\mathrm{log}L)}\displaystyle \frac{{\rm{\Delta }}L}{L}.\end{eqnarray} \tag{ 11 }$$

The total local event rate density is therefore

$$\begin{eqnarray}{\rho }_{0} & = & {\rho }_{0,\gt {L}_{{\rm{min}}}}={\displaystyle \int }_{{L}_{{\rm{min}}}}^{{L}_{{\rm{max}}}}{\rho }_{0,L}{dL}\\ & \simeq & \displaystyle \sum _{\mathrm{log}{L}_{{\rm{min}}}}^{\mathrm{log}{L}_{{\rm{max}}}}\displaystyle \frac{4\pi }{{\rm{\Omega }}T}\displaystyle \frac{1}{\mathrm{ln}10}\displaystyle \frac{1}{g(L)}\displaystyle \frac{{\rm{\Delta }}N}{{\rm{\Delta }}(\mathrm{log}L)}\displaystyle \frac{{\rm{\Delta }}L}{L},\end{eqnarray} \tag{ 12 }$$

which depends on L_min. Observationally L_min is not well constrained, and one usually adopts the observed minimum luminosity, which is the upper limit of the true L_min. As a result, the derived ρ₀ is in principle only the lower limit of the true value. To be specific, throughout the paper, we always specify a minimum luminosity whenever an event rate density is quoted.

With the definition of ρ₀ (Equation (12)), the specific event rate density can also be written as

$$\begin{eqnarray}&&{\rho }_{0,L}={\rho }_{0}{\rm{\Phi }}(L).\end{eqnarray} \tag{ 13 }$$

Within the framework that the LF does not depend on redshift (the approach adopted in this paper), the redshift-dependent event rate density can be written as

$$\begin{eqnarray}&&\rho (z)={\rho }_{0}f(z),\end{eqnarray} \tag{ 14 }$$

where f(z) is the redshift evolution function, the form of which depends on the properties of the transients.

The LF of a certain type of transient can be derived by displaying the specific event rate density ρ_0,L ∝ Φ(L) as a function of L. By separating the data into different luminosity bins, we use the observed numbers to map the relevant ρ_0,L and then fit the data points by several empirical model forms. The simplest model is a single-power-law (SPL) form, i.e.,

$$\begin{eqnarray}&&{\rm{\Phi }}(L)\propto {L}^{-\alpha }.\end{eqnarray} \tag{ 15 }$$

If the model does not define the data well, we introduce a smoothly connected BPL form

$$\begin{eqnarray}&&{\rm{\Phi }}(L)\propto {\left[{\left(\displaystyle \frac{L}{{L}_{{\rm{b}}}}\right)}^{\omega {\alpha }_{1}}+{\left(\displaystyle \frac{L}{{L}_{{\rm{b}}}}\right)}^{\omega {\alpha }_{2}}\right]}^{-1/\omega },\end{eqnarray} \tag{ 16 }$$

where α₁ and α₂ are the PL indices before and after the break luminosity L_b, and ω defines the sharpness of the break. In more complicated cases (e.g., the GLF of LGRBs), one needs another PL segment to fit the data, and we introduce a triple-power-law (TPL) form, i.e.,

$$\begin{eqnarray}{\rm{\Phi }}(L) & \propto & \left\{{\left[{\left({\left(\displaystyle \frac{{L}_{{\rm{b}},1}}{{L}_{{\rm{b}},2}}\right)}^{{\omega }_{1}{\alpha }_{2}}+{\left(\displaystyle \frac{{L}_{{\rm{b}},1}}{{L}_{{\rm{b}},2}}\right)}^{{\omega }_{1}{\alpha }_{3}}\right)}^{-1/{\omega }_{1}}{\left(\displaystyle \frac{L}{{L}_{{\rm{b}},1}}\right)}^{-{\alpha }_{1}}\right]}^{-{\omega }_{2}}\right.\\ & & +{\left.{\left[{\left({\left(\displaystyle \frac{L}{{L}_{{\rm{b}},2}}\right)}^{{\omega }_{1}{\alpha }_{2}}+{\left(\displaystyle \frac{L}{{L}_{{\rm{b}},2}}\right)}^{{\omega }_{1}{\alpha }_{3}}\right)}^{-1/{\omega }_{1}}\right]}^{-{\omega }_{2}}\right\}}^{-1/{\omega }_{2}},\end{eqnarray} \tag{ 17 }$$

where α₁, α₂, and α₃ are the PL indices for three segments, L_b,1 and L_b,2 are the two break luminosities, and ω₁ and ω₂ are the sharpness parameters at the two breaks.

GRBs have a large enough sample to study the redshift evolution effect. We investigate the evolution effect for LGRBs and SGRBs separately. For each class, we first separate the observed GRBs into several redshift bins and then apply our method to map the corresponding local specific event rate density using the GRBs in that specific redshift bin only.⁵ Practically, for a redshift bin (z₁, z₂) around a certain redshift z, we change the integration limits in Equation (2) to z₁ and min(z₂, z_max(L)), respectively, and repeat the procedure laid out in Section 2.1. The derived local specific event rate density and event rate density are expressed as ${\rho }_{0,L}^{z}$ and ${\rho }_{0,\gt {L}_{m}}^{z},$ respectively, denoting that they are derived in the redshift bin around z. Notice that ${\rho }_{0,L}^{z}$ still denotes the local value. By applying a proper correction with the redshift evolution function f(z) (see details in Section 3), one can obtain the LF in the redshift bin (z₁ < z < z₂). If GRB LF does not evolve with redshift, then the results derived from different z bins should remain the same. To minimize the truncation effect by the flux limits of the detectors, we also use a subsample with a higher threshold (with peak photon flux larger than 1.8 photons s⁻¹ cm⁻²) in the derivations and compare the results with two thresholds.

Alternatively, we also repeat the analysis under the same assumption adopted in some recent papers (e.g., Pescalli et al. 2015; Petrosian et al. 2015; Yu et al. 2015), i.e., fix the shape of the LF and assume that L_b evolves as a PL with k ∼ 2.3. We can then map the LF from the data by correcting the luminosity L at z to the "local" value, i.e., ${L}_{0}=L/{(1+z)}^{2.3}.$ Applying the method in Section 2.1 using L₀ instead of L would lead to the "local" LF.

To perform our analysis, the redshift of an event is needed. For GRBs, not all events have redshift measurements. In order to properly account for their event rate density, one needs to correct the derived values based on the z-known sample by the ratio between the total number and the z-known event number. For HL-LGRBs and SGRBs, this correction factor is approximately 3, and we adopt it in our derivations.

The detected rates also depend on the detector's spectral window. For GRBs, BATSE observations suggested that the short-to-long ratio is about 1:3 (Paciesas et al. 1999), whereas Swift BAT, which has a softer bandpass, only gets a 1:10 ratio (Sakamoto et al. 2011; Qin et al. 2013). Since Swift can essentially detect all BATSE LGRBs, one needs to correct for the Swift-detected SGRBs by another factor of ∼3 for the inferred SGRB event rate density. This factor has been taken into account in our derivations.

For the other types of transients, the samples are not large enough to access whether we have missed some events owing to the imperfect instrumental spectral window or the lack of redshift measurements. We therefore derive the event rate density using the observed events only.

The three important instrumental parameters that are relevant to our derivations are the operation time (T), the field of view (Ω), and the flux sensitivity threshold (F_th). Table 1 lists the three parameters of the high-energy detectors used to study various high-energy transients discussed in this paper. While the first two parameters are straightforwardly defined, the definitions of sensitivity thresholds are nontrivial. This is particularly true for wide-field triggering detectors such as Swift BAT, with which most of the GRBs and jetted TDEs were detected (e.g., Lien et al. 2014 for a detailed description of the Swift BAT trigger algorithm). In this paper, we adopt an approximate threshold for each detector. For Swift-detected events, we also adopt slightly different values for different types of events.

Most GRBs were detected through the rate trigger algorithm by BAT, whereas some low-flux events, such as the LL-LGRB 060218, were detected through an image trigger algorithm. The image trigger is an additional trigger algorithm to accumulate photons from a source in a comparably longer time to look for transients that are not bright enough to make a rate trigger.

For the rate-trigger GRB events detected by Swift/BAT, we adopt slightly different values for different subtypes. In order to have the redshift of a GRB measured, the burst usually needs to have a bright enough optical afterglow. On average, the z-known GRBs are brighter and therefore have a higher flux threshold than the standard rate trigger flux threshold. Based on the lowest value of the flux distributions of our subsample, we adopt F_th = 3 × 10⁻⁸ erg cm⁻² s⁻¹ for rate-triggered LGRBs (including HL-LGRBs and LL-LGRB 080517) and 10⁻⁷ erg cm⁻² s⁻¹ for SGRBs.

Taking HL-LGRBs as an example, we derive the flux threshold based on the lower end of the observed photon flux distribution. We adopt the photon flux 0.3 photons cm⁻² s⁻¹ as the threshold.⁶ The transformation from photon flux to specific flux (Equation (26)) requires the information of the spectrum. We consider a typical Band function spectrum (Band et al. 1993) with α = −1, β = −2.3, E_peak = 200 keV at redshift z = 1. This threshold photon flux is translated to F_th = 3 × 10⁻⁸ erg cm⁻² s⁻¹.

The LL-LGRBs 060218 and 100316D and the two Swift-detected TDE jets were detected through image trigger. In the case of image trigger algorithm, the threshold flux depends on the trigger duration T_td, with the dependence ${F}_{{\rm{th}}}\propto {T}_{{\rm{td}}}^{-1/2}.$ The trigger duration of GRB 060218 was about 80 s,⁷ with a mean flux of ∼2.8 × 10⁻⁹ erg cm⁻² s⁻¹ (T. Sakamoto 2015, private communication). This roughly corresponds to a threshold flux

$$\begin{eqnarray}&&{F}_{{\rm{th}}}=2.5\times {10}^{-8}{T}_{{\rm{td}}}^{-1/2}\;\mathrm{erg}\;{\mathrm{cm}}^{-2}\;{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 18 }$$

which we adopt to also calculate the threshold flux for other image trigger events. The trigger duration of GRB 100316D is 64 s.⁸ The two jetted TDE events Sw J1644+57 and Sw J2058+05 had a trigger duration of 64 s⁹ (Burrows et al. 2011) and 4 days (Cenko et al. 2012), respectively, and the corresponding F_th are used to derive ρ_0,L of jetted TDEs. In fact, Sw J1644+57 was image-triggered four times. The 64 s trigger duration was the relevant one at the peak luminosity.

XRO 080109 was serendipitously discovered by Swift XRT. We use a count rate of 0.03 counts s⁻¹ for XRT, with which source variability can be detected for a ∼1000 s observation. This corresponds to a flux threshold F_th ∼ 10⁻¹² erg cm⁻² s⁻¹.

Normal X-ray TDEs were usually discovered from archival survey data of various X-ray missions, such as the ROSAT PSPC All-Sky Survey (RASS; Voges et al. 1999), the XMM-Newton Slew Survey Source Catalog (XMMSL1; Saxton et al. 2008), and the Chandra ACIS archival data. The exposure-time-dependent flux sensitivity thresholds of these three detectors are listed in Table 1. For each TDE event, we consider the real exposure time to determine its F_th.

LGRBs (both HL and LL) and SBOs are associated with the deaths of massive stars. To first order, their redshift distribution should track the history of star formation. There is evidence that at high z, the GRB rate may exceed what the star formation history (SFH) predicts (e.g., Kistler et al. 2008; Li 2008; Qin et al. 2010; Virgili et al. 2011a; Robertson & Ellis 2012). However, for a wide redshift span, the SFH is a good proxy of the redshift distribution of LGRBs. In this paper, we adopt the rough analytical model of SFH derived by Yüksel et al. (2008) using the observational data:

$$\begin{eqnarray}{f}_{{\rm{LGRB/SBO}}}(z) & = & \left[{(1+z)}^{3.4\eta }+{\left(\displaystyle \frac{1+z}{5000}\right)}^{-0.3\eta }\right.\\ & & +{\left.{\left(\displaystyle \frac{1+z}{9}\right)}^{-3.5\eta }\right]}^{\frac{1}{\eta }},\end{eqnarray} \tag{ 19 }$$

where η = −10. At z < 4, this function is directly derived from the SFH inferred from the UV and far-IR galaxy data (Hopkins & Beacom 2006), which is independent of the GRB observations. At z = 5–7, the SFH is enhanced from the galaxy-constrained SFH to compensate the observed GRB excess at high z, which can be explained by the deficiency of the observed low-luminosity star-forming galaxies missed in surveys but traced by GRBs. We adopt Equation (19) to study both LGRBs and SBOs. The latter are only observed in the nearby universe, so that the modification at high z does not enter the problem.

Unlike LGRBs, most SGRBs do not directly trace SFH. Observations suggest that most SGRBs are consistent with having an origin not related to massive star deaths. The leading scenario is mergers of double compact star systems, e.g., two neutron stars (NS–NS) or a neutron star and a black hole (NS–BH; Gehrels et al. 2004; Barthelmy et al. 2005; Fox et al. 2005; Fong et al. 2010; see Berger 2014 for a recent review).

In order for a merger to occur, a compact star binary system needs to go through a long inspiral phase defined by energy loss of the system through gravitational wave radiation. The redshift distribution of SGRBs therefore needs to account for an additional time delay due to inspiral with respect to the creation of the compact binary system, which itself traces the SFH. The distribution of the merger delay timescale is unfortunately not known. Practically, one assumes some empirical forms of the merger delay timescale distribution models and applies the data to derive best parameters for the delay model. Three types of merger delay time models have been discussed in the literature (e.g., Virgili et al. 2011b; Wanderman & Piran 2014 and references therein): PL delay model, Gaussian delay model, and lognormal delay model. Current data support either a Gaussian (Virgili et al. 2011b) or lognormal (Wanderman & Piran 2014) delay model, with the (PL) model disfavored (even though not completely ruled out; Virgili et al. 2011b; Wanderman & Piran 2014). Table 2 lists the three models with the best parameters currently constrained by the SGRB data.

Table 2.  Best-fit Merger Delay Models of SGRBs with Respect to Star Formation History

Delay Model	Formula	Best-fit Parameters	Reference
Gaussian (G)	${m}_{{\rm{G}}}(\tau )d\tau =\mathrm{exp}\left(-\frac{{(\tau -{t}_{{\rm{d,G}}})}^{2}}{2{\sigma }_{{\rm{t,G}}}^{2}}\right)/\sqrt{2\pi }{\sigma }_{{\rm{t,G}}}d\tau$	tt,G = 2 Gyr, σt,G = 0.3	(1)
Lognormal (LN)	${m}_{{\rm{LN}}}(\tau )d\mathrm{ln}\tau =\mathrm{exp}\left(-\frac{{(\mathrm{ln}\tau -\mathrm{ln}{t}_{{\rm{d,LN}}})}^{2}}{2{\sigma }_{{\rm{t,LN}}}^{2}}\right)/\left(\sqrt{2\pi }{\sigma }_{{\rm{t,LN}}}\right)d\mathrm{ln}\tau$	tt,LN = 2.9 Gyr, σt,LN = 0.2	(2)
Power law (PL)	${g}_{{\rm{PL}}}(\tau )d\tau ={\tau }^{-{\alpha }_{{\rm{t}}}}d\tau$	αt = 0.81	(2)

References. (1) Virgili et al. 2011b; (2) Wanderman & Piran 2014.

Download table as:  ASCIITypeset image

With the consideration of the merger delay time distribution, it is difficult to construct an analytical model for the redshift distribution of SGRBs. Instead, we perform a series of Monte Carlo simulations based on the SFH and merger delay distribution models to construct several redshift distributions that correspond to the three delay models with the best-fit parameters. First, we randomly generate 10,000 compact star binary systems with a redshift distribution tracking the SFH following the model of Yüksel et al. (2008). Next, we randomly generate the merger delay timescales of all these systems based on the three merger delay timescale models listed in Table 2. For each model, we derive the lookback time of SGRBs by subtracting the merger delay time from the formation time, and we transfer the lookback time to redshift. We repeat the process 10,000 times (each with 10,000 events simulated). By averaging the results, we are able to derive the average redshift distribution of the simulated samples. We fit the derived redshift distribution (for all three merger delay models) using multiple-PL functions and derive an empirical expression of f(z) for each model. The simulated results with best-fit empirical models are shown in Figure 1. The distributions are normalized to unity at the local universe (z = 0). The empirical formulae of f(z) for the three merger delay models are as follows.

Zoom In Zoom Out Reset image size

For the Gaussian delay model (Virgili et al. 2011b), one has

$$\begin{eqnarray}{f}_{{\rm{SGRB}}}^{{\rm{G}}}(z) & = & \left[{(1+z)}^{5.0\eta }+{\left(\displaystyle \frac{1+z}{0.17}\right)}^{0.87\eta }\right.\\ & & +{\left.{\left(\displaystyle \frac{1+z}{4.12}\right)}^{-8.0\eta }+{\left(\displaystyle \frac{1+z}{4.05}\right)}^{-20.5\eta }\right]}^{\frac{1}{\eta }},\end{eqnarray} \tag{ 20 }$$

with η = −2, which is roughly a BPL with redshift breaks at z₁ = 0.45, z₂ = 2.0, and z₃ = 3.0.

For the lognormal delay model (Wanderman & Piran 2014), one has

$$\begin{eqnarray}{f}_{{\rm{SGRB}}}^{{\rm{LN}}}(z) & = & \left[{(1+z)}^{5.7\eta }+{\left(\displaystyle \frac{1+z}{0.36}\right)}^{1.3\eta }\right.\\ & & +{\left.{\left(\displaystyle \frac{1+z}{3.3}\right)}^{-9.5\eta }+{\left(\displaystyle \frac{1+z}{3.3}\right)}^{-24.5\eta }\right]}^{\frac{1}{\eta }},\end{eqnarray} \tag{ 21 }$$

with η = −2, which is roughly a BPL with redshift breaks at z₁ = 0.35, z₂ = 1.5, and z₃ = 2.3.

For the PL model (Wanderman & Piran 2014), one has

$$\begin{eqnarray}{f}_{{\rm{SGRB}}}^{{\rm{PL}}}(z) & = & \left[{(1+z)}^{1.9\eta }+{\left(\displaystyle \frac{1+z}{2.5}\right)}^{-1.2\eta }\right.\\ & & +{\left.{\left(\displaystyle \frac{1+z}{3.8}\right)}^{-4.4\eta }+{\left(\displaystyle \frac{1+z}{7.7}\right)}^{-11\eta }\right]}^{\frac{1}{\eta }},\end{eqnarray} \tag{ 22 }$$

with η = −2.6. This model has a wider redshift distribution compared to the first two models (owing to the wide range of the merger delay time). It is roughly a BPL with redshift breaks at z₁ = 0.42, z₂ = 3.4, and z₃ = 11.3. The SGRB data do not favor this model (Virgili et al. 2011b), even though it is not completely ruled out.

The event rate density of TDEs depends on both the number density of SMBHs and the event rate of TDEs per galaxy. Since the TDE rate of a particular SMBH only depends on the properties of the galaxy itself (e.g., stellar density near the SMBH and the mass of SMBH), on average, it may be reasonable to assume that there is no redshift evolution of the event rate per galaxy. As a result, the f_TDE(z) parameter of TDEs is mostly determined by the evolution of the number density of SMBHs as a function of redshift (e.g., Donnarumma & Rossi 2015). Shankar et al. (2013) constructed the mass density distribution models for SMBHs and active galactic nuclei (AGNs) by considering their growth rate and radiation efficiency. We apply their model to derive f_TDE(z) for TDEs. TDEs can happen only when the tidal disruption radius is larger than the event horizon of SMBHs, which gives an upper limit of the mass of SMBHs for TDEs:

$$\begin{eqnarray}&&{M}_{{\rm{BH}}}\leqslant 1.6\times {10}^{8}\left(\displaystyle \frac{{M}_{*}}{{M}_{\odot }}\right){\left(\displaystyle \frac{{R}_{*}}{{R}_{\odot }}\right)}^{3/2},\end{eqnarray} \tag{ 23 }$$

where M_* and R_* are the mass and radius of the star that is disrupted by the SMBH, both normalized to the solar values. We therefore exclude SMBHs with mass exceeding 10⁸ M_⊙. In the left panel of Figure 2, we present the numerical fits to the number density redshift evolution of SMBHs with two mass ranges (10⁶–10⁷ M_⊙ and 10⁷–10⁸ M_⊙). Assuming a contant TDE rate per galaxy, in the right panel of Figure 2, we present the normalized TDE redshift distribution f_TDE(z) derived from the numetical data based on the model of Shankar et al. (2013). The best-fit empirical model reads

$$\begin{eqnarray}&&{f}_{{\rm{TDE}}}(z)={\left[{(1+z)}^{0.2\eta }+{\left(\displaystyle \frac{1+z}{1.43}\right)}^{-3.2\eta }+{\left(\displaystyle \frac{1+z}{2.66}\right)}^{-7.0\eta }\right]}^{\frac{1}{\eta }},\end{eqnarray} \tag{ 24 }$$

with η = −2. One can see that f_TDE(z) continues to decrease with redshift, reaching 2/5 at z ∼ 1 and ∼10⁻³ at z ∼ 6.

Zoom In Zoom Out Reset image size

Our HL-GRB sample is only limited to Swift GRBs. This is because it is a uniform sample whose size is large enough to derive a well-constrained GLF. We collect all the z-known Swift GRBs before 2014 May 6 (250 HL-LGRBs and 20 SGRBs). This sample consists of more GRBs than previous work by introducing a lower flux threshold, which allow us to better study the features near the low-luminosity end. The data are downloaded from the Swift archival table available at http://swift.gsfc.nasa.gov/archive/grb_table/ (Sakamoto et al. 2008, 2011). For all the bursts, the 1 s peak photon flux and photon index are recorded. For HL-LGRBs, the 1 s peak photon flux is directly derived from their 1 s peak luminosity. For SGRBs, since their durations are typically shorter than 1 s, we apply a photon count rate with a 64 ms resolution to derive the 64 ms peak luminosity. We calculate their 64 ms peak photon flux based on the ratio of the two peak count rates with different temporal resolutions (C_{p,64 ms} and C_{p,1 s}), i.e., ${P}_{64}={P}_{1}({C}_{\text{p,64 ms}}/{C}_{\text{p,1 s}}).$ The 64 ms light curves are from the Swift Burst Ground-analysis Information page (http://gcn.gsfc.nasa.gov/swift_gnd_ana.html), and the 1 s light curves are obtained through rebinning. To ensure the correct match at the peak, the regrouping is such that the time interval at the peak time (64 ms resolution) matches the one for the 1 s peak photon flux provided in the GCN Circular archive. Only a handful of LL-LGRBs were detected so far. Table 3 collects the information of six LL-LGRBs studied in this paper, which were triggered not only by Swift (GRB 060218, GRB 080517, and GRB 100316D) but also by other instruments as well: GRB 980425 by CRGO/BATSE, XRF 020903 by HETE-II,¹⁰ and GRB 031203 by INTEGRAL. The peak photon fluxes of pre-Swift LL-LGRBs are adopted from GCN archives or Sakamoto et al. (2004). The data of Swift LL-LGRBs are also taken from the Swift table.

The time-integrated spectral information is taken from the literature (references provided in Table 3), described by either an SPL with photon index ${\rm{\Gamma }}(N(E)\propto {E}^{-{\rm{\Gamma }}}),$ a Band function characterized by peak energy E_peak and two photon spectral indices α and β (Band et al. 1993), or a PL function with an exponential cutoff (CPL) fit, i.e., $N(E)\propto {E}^{\alpha }\mathrm{exp}(-E/{E}_{c}).$ For the latter two models, an E_peak can be derived from the peak in the νF_ν spectrum. For SPL fits to most BAT spectra (owing to the narrowness of the BAT band), it is believed that the intrinsic spectrum still has a peak energy. With BAT GRBs jointly detected by other wide-band detectors such as Konus/Wind and Fermi/GBM, it was found that there exists a rough correlation between the BAT-band photon index Γ and E_peak, if E_peak is not much beyond the BAT energy band (Zhang et al. 2007; Sakamoto et al. 2009; Virgili et al. 2012). The latest best fit reads (Virgili et al. 2012)

$$\begin{eqnarray}&&\mathrm{log}({E}_{{\rm{peak}}})=(4.34\pm 0.475)-(1.32\pm 0.129){{\rm{\Gamma }}}^{{\rm{BAT}}}\end{eqnarray} \tag{ 25 }$$

with a large scatter, where Γ^BAT is photon index (positive value) defined in the BAT band. We apply this scaling to estimate E_peak for those GRBs whose E_peak is not directly measured. For a consistency check, we have also adopted E_peak values derived by Butler et al. (2007) for a subsample of GRBs (for which E_peak is available from that method). By repeating the calculations, we found that the derived LF using the Butler et al. (2007) method is similar to the LF derived using our method.

For a GRB with peak photon flux P_p, the peak flux can be calculated through

$$\begin{eqnarray}&&{F}_{{\rm{p}}}=\displaystyle \frac{{P}_{{\rm{p}}}{\displaystyle \int }_{{e}_{1}}^{{e}_{2}}{EN}(E){dE}}{{\displaystyle \int }_{{e}_{1}}^{{e}_{2}}N(E){dE}},\end{eqnarray} \tag{ 26 }$$

where N(E) is the photon spectrum of a GRB, which is in the form of the standard Band function (Band et al. 1993)

$$\begin{eqnarray}&&N(E)=A\left\{\begin{array}{l}{\left(\displaystyle \frac{E}{\text{100 keV}}\right)}^{\alpha }\mathrm{exp}(-\displaystyle \frac{E}{{E}_{0}}),\\ \quad E\lt (\alpha -\beta ){E}_{0},\\ {\left[\displaystyle \frac{(\alpha -\beta ){E}_{0}}{\text{100 keV}}\right]}^{\alpha -\beta }\mathrm{exp}(\beta -\alpha ){\left(\displaystyle \frac{E}{\text{100 keV}}\right)}^{\beta },\\ \quad E\geqslant (\alpha -\beta ){E}_{0}.\end{array}\right.\end{eqnarray} \tag{ 27 }$$

Here the integration limits (e₁, e₂) enclose the detector's energy window (e.g., 15–150 keV for Swift BAT). For SGRBs, we use the E_peak data of 13 SGRBs derived by Lü et al. (2015). For the other seven SGRBs whose E_peak was not directly measured, we estimate E_peak using Equation (25). For those GRBs whose Band function parameters are not directly measured, we adopt typical values as α = −1 and β = −2.3 for LGRBs and α = −0.5 and β = −2.3 for SGRBs.

In order to derive the bolometric luminosity (1–10⁴ keV in the cosmological rest frame) from the observed peak flux, we perform a k-correction

$$\begin{eqnarray}&&{L}_{{\rm{p}},\mathrm{bol}}=4\pi {D}_{{\rm{L}}}^{2}{F}_{{\rm{p}}}\cdot k,\end{eqnarray} \tag{ 28 }$$

where D_L is the luminosity distance. The k-correction parameter can be expressed as

$$\begin{eqnarray}&&k=\displaystyle \frac{{\displaystyle \int }_{1/(1+z)}^{{10}^{4}/(1+z)}{EN}(E){dE}}{{\displaystyle \int }_{{e}_{1}}^{{e}_{2}}{EN}(E){dE}}.\end{eqnarray} \tag{ 29 }$$

Massive stars end their lives in catastrophic core collapses when they run out of fuel in the center (Woosley & Weaver 1986). As a massive star undergos core collapse, an outgoing shock surges through the star. When the optical depth of photons trapped in the shock becomes unity, a SBO occurs, which provides the first electromagnetic emission from a supernova event. Before the shock breaks out of the star, only neutrinos and gravitational waves can escape. At the moment of breakout, a short, bright flash is expected, which peaks in ultravoilet or X-rays depending on how compact the star is (Colgate 1975; Klein & Chevalier 1978; Nakar & Sari 2010). The SBO signal therefore carries direct information about the very early stage of core collapse and provides direct constraints on the type of progenitor. Since there is no electromagnetic precursor to alert such an event, detecting an SBO is challenging. In the X-ray and soft γ-ray regime for which our paper focuses on, there are only two confirmed SBOs detected so far. One is the GRB 060218/SN 2006aj association system, which shows an X-ray thermal component with a temperature of ∼0.17 keV in a very long duration (T₉₀ = 2100 s), soft GRB with a smooth light curve (Campana et al. 2006). The other is the XRO 080109/SN 2008D association system, which was serendipitously detected by Swift/XRT on 2008 January 9 (Soderberg et al. 2008). Since no γ-ray counterpart was detected even though this outburst was in the field of view of BAT before and during the burst, a GRB connection was ruled out. These two observations have offered a great opportunity to study the detailed properties about the progenitors. The fact that several other LL-LGRBs seem to share similar properties to GRB 060218 makes some authors suggest that all LL-LGRBs may be associated with SBOs (e.g., Wang et al. 2007; Nakar & Sari 2012).

We use the two confirmed SBO events (shown in Table 4, with very different luminosities) to estimate their event rate densities. The data of GRB 060218/SN 2006aj are already included in the table of LL-LGRBs. The data of XRO 080109/SN 2008D are collected from Soderberg et al. (2008). Its X-ray spectrum was fit by a PL model with a photon index of 2.61 in the band of 0.3–10 keV, which was used for k-correction.

Stellar tidal disruption by an SMBH has been theoretically predicted by Rees (1988). When a star approaches an SMBH, a TDE occurs if the tidal force becomes larger than the self-gravity of the star and if the radius when this happens is outside the black hole event horizon. Part of the disrupted debris falls into the black hole from an accretion disk, giving rise to a bright flare in UV or X-ray band, which lasts from several months to 1 yr. The first TDE was discovered in quiescent galaxy NGC 5905 during the RASS (Bade et al. 1996; Komossa & Bade 1999), which showed the characteristic luminosity decay law L ∝ t^−5/3 expected for fallback accretion. Later, several more TDEs were detected from RASS, namely, RX J1242, RX J1624, and RX J1420 (Grupe et al. 1999; Komossa & Greiner 1999; Greiner et al. 2000). In the recent decade, more TDE candidates have been discovered by XMM-Newton and Chandra, mostly by comparing survey catalogs with the archival data (Esquej et al. 2007, 2008; Maksym et al. 2010; Saxton et al. 2012; Maksym et al. 2013), as well as by serendipitous detections (Lin et al. 2011). Another ROSAT source was identified as a TDE candidate as it disappeared in the subsequent observations with XMM-Newton and Chandra (Cappelluti et al. 2009). Right now about a dozen X-ray TDE candidates have been discovered. All these TDEs have large amplitudes and soft X-ray spectra, whose host galaxies show no sign of AGN activity (Komossa 2012). Their observed maximum luminosities range from 10⁴² to 10⁴⁵ erg s⁻¹.

Recently, two special TDEs were detected by Swift. These two TDEs, i.e., Swift J1644+57 (Bloom et al. 2011; Burrows et al. 2011) and Swift J2058.4+0516 (Cenko et al. 2012), showed some distinct features. The peak luminosity of Sw J1644+57 was around 10⁴⁸ erg s⁻¹, and the event was followed by a radio counterpart (Zauderer et al. 2011). Owing to the super-Eddington nature of the events, these TDEs have been interpreted as relativistic jets launched from the central black hole (Bloom et al. 2011; Burrows et al. 2011). It has been claimed that there is a low probability that normal TDEs also host a jet similar to Swift J1644+57 (van Velzen et al. 2013). However, it remains unknown whether two types of TDEs are indeed intrinsically different from others, and if so, what could be the main reason to make the difference. By measuring the spin parameter of the central black holes of the two Swift TDEs within the theoretical framework of the Blandford–Znajek mechanism, Lei & Zhang (2011) found that both black holes carry a moderately high spin. They then suggested that black hole spin may be the key factor to make the dichotomy of TDEs, and only black holes with rapid spin can launch relativistic jets during TDEs. Further modeling suggests that the jet model can successfully account for the X-ray (Lei et al. 2013; Tchekhovskoy et al. 2014) and the radio (Metzger et al. 2012; Wang et al. 2014; Liu et al. 2015) data.

Table 5 lists the data for all the TDE candidates. TDE flares last much longer than GRBs and SBOs. They usually have a relatively fast rising phase, reach and stay at the peak for some time, and then decay roughly with a PL L ∝ t^−5/3. The luminous state usually lasts from several months to 1 yr (Rees 1988). In our sample, only two events may have been detected in both the rising and declining phases, so that the peak luminosity was measured (Komossa & Bade 1999; Esquej et al. 2008). For other TDEs, one did not detect the sources in both the rising and the declining phases, so that the peak luminosity cannot be well constrained. The maximum luminosity during the detected phase is only the lower limit of the peak luminosity. However, since the light-curve peak is rounded and spreads in several weeks (Lodato et al. 2009), the observed peak luminosity may not be too different from the true peak luminosity.

Owing to the narrow bandpass of the XRTs that detect TDEs, the spectral shape for TDEs is not well constrained. Except for the BAT-detected TDEs for which a treatment similar to GRBs can be applied, for the majority of TDEs, we only apply an empirical relation to estimate the bolometric luminosity by multiplying X-ray luminosity by a factor of 1.4, i.e. (Maksym et al. 2010),

$$\begin{eqnarray}&&{L}_{{\rm{p}},\mathrm{bol}}=1.4\times {L}_{{\rm{p}}}.\end{eqnarray} \tag{ 30 }$$

For Swift TDEs, a PL spectrum in the BAT band is reported for both events with a photon index of 1.8 and 1.6, respectively (Burrows et al. 2011; Cenko et al. 2012). Even though the spectrum of jetted TDEs is not known, we speculate that they have a nonthermal spectrum with E_p not far above the BAT band. As a result, we apply Equation (25) to estimate E_p and apply Equation (29) to estimate the k-correction factor assuming the standard Band function parameters for the spectrum. This led to a k-correction factor of 2.1 and 2.2 for Sw J1644+57 and Sw J2058+05, respectively.

Since LGRBs and SGRBs have different physical origins (massive star core collapse versus compact star mergers), we derive their event rate density and LF separately. For each type, we first derive the GLF by ignoring the possible redshift evolution effect. Then we dedicate one subsection to discussing the possible evolution effect. Within LGRBs, the LL-LGRBs have been claimed to form a distinct component in the LF (Liang et al. 2007; Virgili et al. 2009), which may have a somewhat different physical origin (e.g., Liang et al. 2007; Bromberg et al. 2011; Nakar & Sari 2012). In our analysis, we adopt two approaches. First, we derive the LF of HL-LGRBs and LL-LGRBs separately (e.g., Liang et al. 2007; Virgili et al. 2009). In the second approach, we derive the LF of LGRBs jointly by fitting the data together with a two-component (i.e., TPL) LF. This approach is justified in view of the fact that both LL-LGRBs and HL-LGRBs are associated with SNe Ic and therefore may share a common physical origin.

5.1.1. Global LF

Figure 3 shows the results of GRBs. The event rate density above a particular L as a function of the bolometric luminosity L is presented in the left panel, and the specific event rate density as a function of L, which describes the GLF, is presented in the right panel. In both panels, the HL-LGRBs (red), LL-LGRBs (blue), and SGRBs (black) are presented separately. The luminosity bin is taken as 0.3 (LGRBs) or 0.7 (SGRBs) in the logarithmic space. The horizontal errors denote the luminosity bins, whereas the vertical errors are calculated from small-sample statistics (Gehrels 1986). The best-fit lines for all three subtypes of GRBs are also presented. The fitting results are summarized in Tables 6 and 7. For all the fitting parameters, we also present the 1σ range of errors based on 5000 sets of Monte Carlo simulations.

Zoom In Zoom Out Reset image size

Table 6.  The Best-fit Luminoity Functions of Different Types of Extragalactic High-energy Transients

Type	Fit Model	α1	α2	α3	Lb (erg s−1)
HL-LGRBs	TPL	${2.2}_{-0.2}^{+0.4}$	${1.0}_{-0.1}^{+0.1}$	${2.0}_{-0.3}^{+0.3}$	${5.0}_{-1.3}^{+3.0}\times {10}^{50},$ ${7.1}_{-3.0}^{+4.3}\times {10}^{52}$
LL-LGRBs	SPL	${2.3}_{-0.2}^{+0.2}$	...	...	...
Joint HL-/LL-LGRBs	TPL	${1.7}_{-0.1}^{+0.1}$	${1.0}_{-0.1}^{+0.2}$	${2.0}_{-0.3}^{+0.3}$	${1.0}_{-0.3}^{+0.2}\times {10}^{51},$ ${7.8}_{-3.1}^{+2.3}\times {10}^{52}$
SGRBs (G)	SPL	${1.7}_{-0.08}^{+0.08}$	...	...	...
SGRBs (LN)	SPL	${1.6}_{-0.08}^{+0.08}$	...	...	...
SGRBs (PL)	SPL	${1.5}_{-0.08}^{+0.08}$	...	...	...
Joint SBO/LL-LGRB	SPL	${2.0}_{-0.09}^{+0.09}$	...	...	...
TDEs	SPL	${2.0}_{-0.05}^{+0.05}$	...	...	...

Notes. For SGRBs, the results for three merger delay models (Gaussian (G), lognormal (LN), and power law (PL)) are given. The 1σ errors of all the fitting parameters are presented based on 5000 sets of Monte Carlo simulations.

Download table as:  ASCIITypeset image

Table 7.  The Event Rate Density of Various Transients Given an Observed Minimum Luminosity Threshold and a Typical Luminosity Threshold

Type	Lm (erg s−1)	${\rho }_{0,\gt {L}_{m}}({\mathrm{Gpc}}^{-3}\;{\mathrm{yr}}^{-1})$	L' (erg s−1)	${\rho }_{0,\gt {L}^{\prime }}({\mathrm{Gpc}}^{-3}\;{\mathrm{yr}}^{-1})$
HL-LGRBs	3 × 1049	${2.4}_{-0.3}^{+0.3}$	1050	${0.8}_{-0.1}^{+0.1}$
LL-LGRBs	5 × 1046	${164}_{-65}^{+98}$	1046	${440}_{-175}^{+264}$
SGRBs (G)	7 × 1049	${4.2}_{-1.0}^{+1.3}$	1050	${1.3}_{-0.3}^{+0.4}$
SGRBs (LN)	7 × 1049	${3.9}_{-0.9}^{+1.2}$	1050	${1.2}_{-0.3}^{+0.4}$
SGRBs (PL)	7 × 1049	${7.1}_{-1.7}^{+2.2}$	1050	${3.3}_{-0.8}^{+1.0}$
SBOs	1044 a	${1.9}_{-1.2}^{+2.4}\times {10}^{4}$	1047	${14}_{-11}^{+32}$
TDEs	1042	${1.0}_{-0.3}^{+0.4}\times {10}^{5}$	1044	${4.8}_{-2.1}^{+3.2}\times {10}^{2}$
Swift TDEs	1048	${0.03}_{-0.02}^{+0.04}$	...	...

Note.

^aSince there are only two confirmed SBOs, two characteristic luminosities are given around the exact luminosities of the two events.

Download table as:  ASCIITypeset image

The local event rate density of LL-LGRBs is ${164}_{-65}^{+98}\;{\mathrm{Gpc}}^{-3}\;{\mathrm{yr}}^{-1}$ with a minimum luminosity 5 × 10⁴⁶ erg s⁻¹, which is roughly consistent with previous results (Liang et al. 2007; Virgili et al. 2009). The longer working period of Swift until now than the ones relevant for the previous two papers makes the event rate density slightly lower than before. However, it is still around two orders of magnitude higher than that of HL-LGRBs, which is ${2.4}_{-0.3}^{+0.3}\;{\mathrm{Gpc}}^{-3}\;{\mathrm{yr}}^{-1}$ above 3 × 10⁴⁹ erg s⁻¹, or ${0.8}_{-0.1}^{+0.1}\;{\mathrm{Gpc}}^{-3}\;{\mathrm{yr}}^{-1}$ above 10⁵⁰ erg s⁻¹ (a typical luminosity threshold adopted before). This local event rate of HL-LGRBs is slightly lower than Liang et al. (2007) and Wanderman & Piran (2010).

From Equation (13), one can see that ρ_0,L is proportional to the LF Φ(L), with ρ₀ defining the normalization. We use either an SPL or a TPL to fit the ρ_0,L distribution for each subtype of GRBs to characterize their GLFs. The fitting results are shown in Table 6. HL-LGRBs can be fit with a TPL with ${\alpha }_{1}^{{\rm{HL}}}={2.2}_{-0.2}^{+0.4},$ ${\alpha }_{2}^{{\rm{HL}}}={1.0}_{-0.1}^{+0.1},$ ${\alpha }_{3}^{{\rm{HL}}}={2.0}_{-0.3}^{+0.3}$ and the break luminosity ${L}_{{\rm{b}},1}^{{\rm{HL}}}={5.0}_{-1.3}^{+3.0}\times {10}^{50}\;\mathrm{erg}\;{{\rm{s}}}^{-1},$ ${L}_{{\rm{b}},2}^{{\rm{HL}}}={7.1}_{-3.0}^{+4.3}\times {10}^{52}\;\mathrm{erg}\;{{\rm{s}}}^{-1}.$ The luminosity distribution of LL-LGRBs can only be fit using an SPL with ${\alpha }^{{\rm{LL}}}={2.3}_{-0.2}^{+0.2}.$ The LL-LGRBs are not the straightforward extension of HL-LGRBs to low luminosities. But the slope of LL-LGRB is similar to that of the first component of HL-LGRBs. The normalization of LL-LGRBs is a little bit lower as we apply a lower threshold for LL-LGRBs than HL-LGRBs. Therefore, it may be possible that LL-LGRBs follow the extension of LF of HL-LGRBs if we take a TPL fit to the joint LL and HL samples. Such a fit is presented in Figure 4, with the best-fit parameters being ${\alpha }_{1}={1.7}_{0.1}^{+0.1},$ ${\alpha }_{2}={1.0}_{-0.1}^{+0.2},$ ${\alpha }_{3}={2.0}_{-0.3}^{+0.3}$ and the break luminosities ${L}_{{\rm{b}},1}={1.0}_{-0.3}^{+0.2}\times {10}^{51}\;\mathrm{erg}\;{{\rm{s}}}^{-1},$ ${L}_{{\rm{b}},2}\;={7.8}_{-3.1}^{+2.3}\times {10}^{52}\;\mathrm{erg}\;{{\rm{s}}}^{-1}.$ This is similar to the TPL fit to HL-LGRBs alone except for a slightly shallower α₁, which is compromised by the slight mismatch between LL-LGRBs and HL-LGRBs.

Zoom In Zoom Out Reset image size

The TPL nature of the joint GLF for LGRBs is interesting. Such a feature was only noticeable in the past when LL-LGRBs are included (Liang et al. 2007; Virgili et al. 2009). With the current sample, we find that it is required even for HL-LGRBs alone. With the joint fit LL- and HL-LGRBs, we find that the steep component in the low-luminosity end now extends all the way to ∼10⁵¹ erg s⁻¹, so that no clear separation between the LL- and HL-LGRB population is seen.

We would like to stress that if we focus on the two high-L segments in the GLF, our results (α₂ ∼ 1.0, α₃ ∼ 2.0, L_b,2 ∼ 7.8 × 10⁵² erg s⁻¹) are broadly consistent with previous results: α^HL ∼ 0.8, β^HL ∼ 2.6, and L_b^HL ∼ 2.5 × 10⁵² by Liang et al. (2007), and α^HL ∼ 1.17, β^HL ∼ 2.44, and ${L}_{{\rm{b}},2}^{{\rm{HL}}}\sim 3.1\;\times {10}^{52}\;\mathrm{erg}\;{{\rm{s}}}^{-{\rm{s}}}$ by Wanderman & Piran (2010).¹¹

According to Equations (1) and (11), the indices of ${\rho }_{0,\gt L}$ and ρ_0,LdL should be roughly the same, i.e., the index of ${\rho }_{0,\gt L}$ should be greater than the index of ρ_0,L by one. This is generally satisfied for most of the transients studied in this paper (see all the indices marked in Figures 3–8).

5.1.2. Luminosity Function Evolution

Using the method laid out in Section 2.2, we study the evolution effect of LGRB LF. The results of ${\rho }_{0,L}^{z}$ are shown in the left panel of Figure 5. For each redshift bin, we fit the LF with an SPL or BPL if the latter is needed. The fitting parameters are presented in Table 8. One can see that indeed there is an apparent LF evolution effect. However, there is no clear pattern to quantify the evolution. The right panel of Figure 5 shows the redshift-dependent break luminosity. In some redshift bins a break is clearly seen. However, in some other bins, the break either does not exist or is simply required by only one data point with low significance. For SPL fits, we place either a lower limit or an upper limit based on the highest- or lowest-luminosity bin. As shown in the right panel of Figure 5, there is no clear pattern to quantify the evolution effect. Furthermore, the PL indices α₁ and α₂ also show significant evolutions (variations) in different redshift bins. Since different redshift bins have different L_m, and since the behavior below L_m is poorly constrained by the data, in Table 8 we choose different L_m for different redshift bins. For the nearby universe, we also get a ${\rho }_{0}^{z}$ at higher L_m. We obtain a consistency of ${\rho }_{0}^{z}$ derived from data from different redshift intervals.

Zoom In Zoom Out Reset image size

Table 8.   Best-fit LF Parameters in Different Redshift Bins for the Full Sample and Subsample (P_p > 1.8 photons cm⁻² s⁻¹)

Full Sample
Redshift	α1	α2	Lb (erg s−1)	${\rho }_{0,\gt {L}_{m}}^{z}$ (Gpc−3 yr−1)a	Lm (erg s−1)
0 < z < 1	1.6	1.8	4.0 × 1051	${1.6}_{-0.2}^{+0.3}$ (${0.3}_{-0.1}^{+0.1}$)	1050 (1051)
0.5 < z < 1.5	1.4	1.7	≥2.0 × 1053, or ≤5.1 × 1049	${1.2}_{-0.2}^{+0.2}$ (${0.4}_{-0.1}^{+0.1}$)	1050 (1051)
1 < z < 3	1.1	1.6	1.2 × 1052	${0.4}_{-0.1}^{+0.1}$	1051
2 < z < 4	1.5	>4.0	≥5.5 × 1053, or ≤2.0 × 1051	${0.3}_{-0.1}^{+0.1}$	1052
z > 3	∼0	1.7	≤1.7 × 1052	${0.2}_{-0.1}^{+0.1}$	1053
Subsample
Redshift	α1	α2	Lb (erg s−1)	${\rho }_{0,\gt {L}_{m}}^{z}$ (Gpc−3 yr−1)	Lm (erg s−1)
0 < z < 1	0.1b	1.8	8.5 × 1050	${1.7}_{-0.2}^{+0.3}$ (${0.7}_{-0.1}^{+0.1}$)	1050 (1051)
0.5 < z < 1.5	1.5	...	≥1053, or ≤8.0 × 1050	${0.7}_{-0.1}^{+0.1}$	1051
1 < z < 3	1.1	1.6	3.2 × 1052	${0.6}_{-0.1}^{+0.1}$	1051
2 < z < 4	1.0	1.5	4.3 × 1052	${0.6}_{-0.1}^{+0.1}$	1052
z > 3	∼0	1.9	7.0 × 1052	${0.4}_{-0.1}^{+0.1}$	1053

Notes.

^aThe local event rate density derived from the data in different redshift bins. The minimum luminosity varies at each bin owing to the limited instrument sensitivity. We also give event rate density at the same L_m for comparison. ^bIn the redshift bin 0 < z < 1, LF could be fit by a TPL. Here α₁ and α₂ denote the latter two components and L_b is the second break luminosity.

Download table as:  ASCIITypeset image

To minimize the truncation effect by the flux limit of detectors, we also use a subsample with a higher threshold (with peak photon flux larger than 1.8 photons s⁻¹ cm⁻²); see Figure 6. The subsample consists of fewer GRBs in some redshift intervals (e.g., $z\in (0.5,1.5)$) so that the LF can be fit by an SPL. In the nearby universe, on the other hand, since LL-LGRBs are dropped owing to this high-threshold criterion, the LF demands a TPL fit.¹² The fitting parameters are also presented in Table 8. The evolution of L_b is now better quantified by ${L}_{{\rm{b}}}\propto {(1+z)}^{3.7}$ (Figure 6, right panel), but a signficant variation of α₁ and α₂ values in different redshift bins remains. We therefore conclude that there is no straightforward way to quantify the evolution effect of the LGRB LF.

Zoom In Zoom Out Reset image size

On the other hand, if we assume the evolution law assumed from the latest papers (e.g., Pescalli et al. 2015; Petrosian et al. 2015; Yu et al. 2015), i.e., the PL indices remain unchanged, and only L_b evolves with k ∼ 2.3, we can map the LF from the data by correcting the luminosity L at z to the "local" value, i.e., ${L}_{0}=L/{(1+z)}^{2.3},$ to derive the "local" LF. The result is shown in Figure 7. This local LF can be fit by a BPL with α₁ = 1.5, α₂ = 2.5, and L_b,0 = 51.6 erg s⁻¹. This is consistent with the results presented in previous papers.

Zoom In Zoom Out Reset image size

5.2.1. Global LF

The local event rate density for SGRBs varies slightly for different merger delay models. For a minimum luminosity 7 × 10⁴⁹ erg s⁻¹, it is ${4.2}_{-1.0}^{+1.3},$ ${3.9}_{-0.9}^{+1.2},$ ${7.1}_{-1.7}^{+2.2}\;{\mathrm{Gpc}}^{-3}\;{\mathrm{yr}}^{-1}$ for the Gaussian, lognormal, and PL delay models, respectively. Taking a typical minimum luminosity 10⁵⁰ erg s⁻¹ as adopted by previous authors, the numbers are ${1.3}_{-0.3}^{+0.4},$ ${1.2}_{-0.3}^{+0.4},$ and ${3.3}_{-0.8}^{+1.0}\;{\mathrm{Gpc}}^{-3}\;{\mathrm{yr}}^{-1},$ respectively. The local event rate density for the lognormal model is slightly lower than the value reported by Wanderman & Piran (2014), as we use a slightly lower flux threshold for Swift BAT-detected SGRBs.

SGRBs come with a much smaller number of events than LGRBs. Assuming that all the SGRBs with redshift measurements are of a compact star merger origin, we derive their LF in Figure 8. The three different merger delay models give slightly different results, but in general all three models are consistent with having an SPL LF with an index of ∼1.6 (for details, see Figure 6). This is different from HL-LGRBs, which require a TPL LF. It is also different from Wanderman & Piran (2014), who claimed a BPL LF.

Zoom In Zoom Out Reset image size

One caveat of our treatment is the assumption that all SGRBs are of a compact star merger (Type I) origin. In the Swift era, Zhang et al. (2009) first suspected that some high-L SGRBs at high redshifts may not be of the Type I origin, but may rather originate from massive star core collapse (Type II). They suggested to apply multiwavelength criteria (instead of using γ-ray duration only) to judge the physical origin of a GRB. Later Virgili et al. (2011b) pointed out that the assumption that all SGRBs are of the compact star merger origin (Type I) is disfavored since one cannot simultaneously account for the z − L two-dimensional distribution and the $\mathrm{log}N-\mathrm{log}P$ distribution of SGRBs. They claimed that a good fraction of SGRBs may be of a Type II origin. Bromberg et al. (2013) recently reached the similar conclusion using a different argument based on the duration distribution of GRBs. Wanderman & Piran (2014) used the criteria of Bromberg et al. (2013) and excluded about 1/3 of SGRBs that they suspect to be of a massive star core collapse origin. This may explain the difference between the results of this paper (SPL LF) and that paper (BPL). We did not exclude any SGRB from our sample for the following reason. We believe that any conclusion about the physical category of an SGRB based on the duration information only (e.g., Bromberg et al. 2013) is not reliable. Rather, one should consider multiwavelength criteria (Zhang et al. 2009), especially the host galaxy type and the afterglow location within the host. Host galaxy studies of SGRBs (Fong et al. 2010; Berger 2014) suggested that the hypothesis that all SGRBs belong to the compact star merger (Type I) category is not in conflict with the data. Indeed, some SGRBs excluded by Wanderman & Piran (2014; defined by them as collapsars) actually have large offsets from host galaxies, fully consistent with being due to a compact star merger origin (e.g., GRB 070724 with offset 5.46 ± 0.14 kpc, and GRB 070809 with offset 33.22 ± 2.71 kpc; Fong et al. 2010; Fong & Berger 2013). Furthermore, a recent study of SGRB emission amplitude parameter (Lü et al. 2014) also suggested that most observed SGRBs are not disguised SGRBs due to the "tip-of-iceberg" effect. Even though we believe that there exists a contamination of Type II GRBs in the SGRB sample, without studying the multiwavelength data of SGRBs in detail, we believe that it is more appropriate not to exclude any SGRB in this study.

5.2.2. LF Evolution

For SGRBs that have much less data, it is more difficult to study the possible evolution effect of LF. We manage to divide the SGRB sample into three redshift bins (z < 0.3, 0.3 < z < 0.6, and z > 0.6) and apply the method in Section 2.2 to derive ${\rho }_{0,\gt {L}_{m}}^{z}$ using the data in different redshift bins. The results are shown in Figure 9. The LFs can be all fit with an SPL. The slope in the first redshift bin (0 < z < 0.3), ${1.9}_{-0.3}^{+0.3},$ is somewhat steeper than those of the other two redshift bins (${1.5}_{-0.2}^{+0.2}$ for 0.3 < z < 0.6, and ${1.4}_{-0.2}^{+0.2}$ for z > 0.6). However, the slopes are consistent with each other within error.

Zoom In Zoom Out Reset image size

The results for these two SBO events are presented in Figure 10. In view of the possible connection between SBOs and LL-LGRBs in general, we present LL-LGRBs in the same plot for comparison.

Zoom In Zoom Out Reset image size

Our results suggest that the event rate density is ${3.1}_{-2.0}^{+4.1}\times {10}^{4}\;{\mathrm{Gpc}}^{-3}\;{\mathrm{yr}}^{-1}$ for XRO 080109/SN 2008D-like SBO events (luminosity ∼6.1 × 10⁴³ erg s⁻¹) and is ${11}_{-9}^{+25}{\mathrm{Gpc}}^{-3}\;{\mathrm{yr}}^{-1}$ for GRB 060218/SN 2006aj-like SBO events (luminosity ∼1.5 × 10⁴⁷ erg s⁻¹). The former implies that the local event rate density of SBO is at least 10⁴ times higher than that of HL-LGRBs. A sensitive large field of view X-ray detector would lead to discovery of a large sample of these events.

For both figures we find that XRO 080109/SN 2008D (green data point at low luminosity) follows the extension of LL-LGRBs (blue line). We also perform a joint fit between SBOs and LL-LGRBs and get an SPL LF with slope ∼2.0, which is similar to that of the slope of LL-LGRBs only (α^LL = 2.3). This lends support to the possible connection between LL-LGRBs and SBOs (e.g., Wang et al. 2007; Nakar & Sari 2012; Barniol Duran et al. 2014).

Figure 11 presents the results for TDEs. As mentioned above, two caveats are that the measured luminosities are only the lower limits of the "peak luminosities," and that some uncertainties are associated with the k-correction parameters. Bearing in mind these caveats, the following conclusions can be drawn. The event rate densities at different luminosity bins have a wide distribution, ranging from over 10⁵ Gpc⁻³ yr⁻¹ at 10⁴² erg s⁻¹ to ∼10² Gpc⁻³ yr⁻¹ at 10⁴⁵ erg s⁻¹. In the luminosity range of (10⁴³–10⁴⁴) erg s⁻¹, the event rate density is in the range of (10³–10⁴) Gpc⁻³ yr⁻¹, which is consistent with both theoretical predictions (Wang & Merritt 2004) and the estimates based on observations (Esquej et al. 2008; Luo et al. 2008; Maksym et al. 2010). The event rate density of the two Swift TDEs is $\sim {0.03}_{-0.02}^{+0.04}\;{\mathrm{Gpc}}^{-3}\;{\mathrm{yr}}^{-1}.$ Similar to GRB 060218, these two events were detected through image triggers. In particular, Sw J1644+57 triggered BAT multiple times (all image triggers; Burrows et al. 2011). We use the brightest peak to define the peak luminosity, and we derive the event rate density based on the trigger information of that epoch.

Zoom In Zoom Out Reset image size

The GLF of TDEs (including both normal TDEs and jetted TDEs) can be roughly described by an SPL with α^TDE = 2.0. With the current sample, the jetted TDEs detected by Swift seem to lie in the extension of a normal TDE LF to high luminosities. The event rate density of Sw J1644+57 shows a flattening feature at the highest luminosity. More data are needed to verify whether jetted TDEs form a distinct component in the global TDE LF.

In Figure 12, we plot all high-energy transients in one figure. The left panel shows the local event rate density above a certain luminosity (${\rho }_{0,\gt L}$) as a function of L, and the right panel presents the specific local event rate density (ρ_0,L), which essentially represents the GLF. All the data points are given a 1σ error in vertical axis. The width of the luminosity bin is shown as the horizontal error bar.

Zoom In Zoom Out Reset image size

Intriguingly, all the transients seem to line up to form a rough SPL distribution. A best fit to the GLFs of all the transients gives a slope of α^global = 1.6. All the events lie within the 3σ confidence bounds of the best fit. The region below this correlation line is likely due to an observational bias and could be filled with new types of transients. With a much lower event rate density, these transients may not have been detected within the time span of the modern high-energy astronomy. The region above the correlation line, on the other hand, is not subject to selection effects and must be intrinsic. The existence of such an upper boundary of high-energy transients is intriguing, which may be rooted from more profound physical reasons.