COMPACT OBJECT COALESCENCE RATE ESTIMATION FROM SHORT GAMMA-RAY BURST OBSERVATIONS

Short-duration gamma-ray bursts (SGRBs) are some of the most powerful explosions detected in the universe, releasing intensive bursts of high-energy gamma rays with a peak duration shorter than 2 s (Kouveliotou et al. 1993). The most commonly accepted explanation of their origin is a system of two compact objects, either two neutron stars (NS–NS) or a neutron star and a black hole (NS–BH) coalescing into a BH (Eichler et al. 1989; Narayan et al. 1992; Piro 2005; Rezzolla et al. 2011). Because of the emission of gravitational waves (GWs) during the latest phases of binary evolution, these objects are one of the primary sources for the next generation of ground-based GW detectors such as Advanced LIGO (Smith 2009) and Advanced Virgo (Acernese et al. 2008). Direct detection of a GW signal from a compact binary coalescing system would allow astronomers to gain valuable information on the astrophysics of compact objects, for example, the NS equation of state (Flanagan & Hinderer 2008; Read et al. 2009), as well as probe fundamental physics by testing the Lorentz invariance principle (Ellis et al. 2006) and general relativity in the strong-field regime (Will 2005), or by setting limits on the graviton mass (Stavridis & Will 2009; Keppel & Ajith 2010). Direct detection of a GW signal coinciding with a gamma-ray burst (GRB) optical counterpart could provide additional insights on astrophysics and even cosmology. The measure of the redshift of a GRB coinciding with a GW detection could allow astronomers to directly determine the distance of the system (see, e.g., Nissanke et al. 2010 and references therein). Coincident detections could also significantly improve the determination of the Hubble parameter by GW observations (Dalal et al. 2006; Del Pozzo 2011).

In this context, it is crucial to have reliable estimates of the local merger rate of compact objects and the number of expected coincident GW-SGRB observations in the advanced GW detector era. In this paper, we present a simple estimate of these quantities by using SGRB data from SWIFT observations.³ In order to avoid selection bias, we calculate the number of expected coincident observations by restricting the sample of SWIFT data to observations with determined redshift and certain association to an optical counterpart. In contrast to a previous study by one of the authors (Dietz 2011), we also include the observed GRB luminosities in the analysis.

The paper is organized as follows. In Section 2, we define the SWIFT data sample and discuss the theoretical model that is used to fit the observations. In Section 3, we present the results and compare them to other published rate estimates. The Appendix contains details on the rate functions used in the analysis.

We restrict our analysis to a set of SGRBs with reliable redshift measurement, i.e., to SGRBs that can be associated with a galaxy of known spectroscopic redshift with a high probability of being the host galaxy of the SGRB. We omit from the analysis SGRBs without an observed optical afterglow and SGRBs that are not associated with a host galaxy within the error circle of the observation. This allows us to remove any instrumental bias with respect to SGRBs detected by other missions. Table 1 shows a list of the 14 SWIFT SGRBs that pass our selection criteria.

The luminosity of these SGRBs can be computed using their redshift and fluence information. The fluence S is divided by the SGRB duration to estimate the flux F in the relevant frame for the detection threshold of the satellite (the observer frame). Since the observed fluence depends on the spectral properties of the source and the energy response of the detector (15–150 keV for the Burst Alert Telescope (BAT) instrument on board SWIFT; see Barthelmy et al. 2005), and the observations are over cosmological distances, two identical sources at different distances may show a spectral shift and a change of fluence. Expressing the observed photon number spectrum with the Band function (Band et al. 1993), this spectral shift can be calculated as a function of the redshift (Cao et al. 2011). Since most SGRB sources have redshift smaller than z = 1, the effect of the spectral shift is less than 10% (Cao et al. 2011). This error is small compared to other statistical and systematical errors and will be neglected in our analysis. The apparent luminosity of the SGRB is

$$\begin{equation} L= 4\pi d_{\mathrm{L}}^{2}(z)\,F \simeq 4\pi d_{\mathrm{L}}^{2}(z)\frac{S}{T_{90}}\,, \end{equation} \tag{ 1 }$$

where d_L(z) is the luminosity distance for a given redshift z, F is the mean flux, S is the measured fluence, and T₉₀ is the time over which the burst emits 90% of its total energy. Throughout this paper, we consider a standard flat-Λ cosmology with H₀ = 71 km s⁻¹, Ω_M = 0.27, and Ω_Λ = 0.73. The left panel of Figure 1 shows the distribution of the observed luminosities as a function of the redshift. The solid line indicates the approximate detector's sensitivity threshold:

$$\begin{equation} d_\mathrm{max}(L) = \sqrt{\frac{L}{4\pi F_\mathrm{lim}}}\,, \end{equation} \tag{ 2 }$$

where the flux threshold is F_lim = 5 × 10⁻⁹ erg s⁻¹ cm⁻² (Cao et al. 2011). The gray area in Figure 1 defines the so-called redshift desert, a region between z ≃ 1 and z ≃ 2 where spectroscopic redshift determinations are difficult to obtain (Fiore et al. 2007). In the following, we choose a conservative approach and further restrict our data sample to z < 1, leaving 13 data points for our analysis.

Zoom In Zoom Out Reset image size

To determine the local merger rate, we fitted the sample with several commonly used luminosity functions (Dietz 2011; Chapman et al. 2009; Guetta & Stella 2008). The lognormal function and the Schechter function were found to provide acceptable fits. The cumulative distribution of the SGRBs as a function of the luminosity and the lognormal and Schechter functions are shown in Figure 1. Since the lognormal fit is slightly better than the Schechter fit, we restrict our analysis to the lognormal luminosity function

$$\begin{equation} \phi (L)\propto \frac{1}{L}\mathrm{exp}\left(-\frac{\mathrm{log}L-\mathrm{log}L_{0}}{2\sigma ^{2}}\right)\,, \end{equation} \tag{ 3 }$$

where L₀ is the mean (peak) value of the luminosity and σ is the width of the distribution. These two parameters are determined by fitting the function to the 13 SGRBs in the sample. For the sake of simplicity, we do not consider evolutionary effects on the luminosity function. Although it is reasonable to assume that these effects are small compared to other statistical and systematical errors, some physical processes depend on the metallicity of the progenitors, which is a function of the redshift (Belczynski et al. 2010, 2011). The right panel of Figure 1 shows the cumulative distribution of SGRBs as a function of luminosity and the lognormal fit in Equation (3).

Since the fit is made on observed SGRBs, the luminosity function ϕ'(L) must be rescaled to the volume where SWIFT is sensitive. Assuming an isotropic distribution of SGRBs, we write

$$\begin{equation} \phi (L) \propto \phi ^{\prime }(L)/ d_{\mathrm{max}}^{3}(L)\,, \end{equation} \tag{ 4 }$$

where $d_{\mathrm{max}}^{3}(L)$ is the maximum luminosity distance where an SGRB can be detected by SWIFT. The number of observable SGRBs within a redshift distance z is

$$\begin{equation} N^{\prime }(z)=N_{0}\int _{0}^{z}dz^{\prime }\frac{R(z^{\prime })}{1+z^{\prime }}\frac{dV(z^{\prime })}{dz^{\prime }}\int _{L_{\mathrm{min}}(z^{\prime })}^{\infty }\phi (L)dL\,, \end{equation} \tag{ 5 }$$

where dV(z')/dz' is the comoving volume element and N₀ is a normalization factor. The rate function R(z) describes the formation rate of the binary systems per comoving volume as a function of the redshift. Since a binary system of compact objects is formed from massive progenitor stars, R(z) may be assumed to follow the star formation rate. The time difference from the formation of the compact objects to the coalescence of the binary is likely on the order of the Gyr (Belczynski et al. 2006). Thus, there is a significant delay with respect to the star formation rate. This is taken into account by using a delayed rate function in Equation (5). Figure 2 shows the different rate functions that are used in our analysis. The Appendix contains explicit expressions and references.

Zoom In Zoom Out Reset image size

Since Equation (5) describes the number of SGRBs that may be potentially observed, N'(z) must be rescaled to fit the number of SGRBs used in the analysis. From 2005 through 2011 (T ≃ 6 yr), SWIFT observed 46 SGRBs. (We neglect any downtime of the satellite due to technical issues or other constraints.) Thus, N'(z) must be divided by the factor f_R = 14/46, with an estimated error of $\sqrt{14}/46\sim 10$%. The SWIFT field of view is about 1.4 sr (Barthelmy et al. 2005), corresponding to a visible fraction of the sky approximately equal to f_FOV ≃ 10%. Observations have shown that about 15% of all SGRBs may be created by soft gamma repeaters (SGRs; Nakar et al. 2006; Chapman et al. 2009). Since SGRs are typically less bright than ordinary SGRBs, events at larger redshifts might be composed mainly of SGRBs. If we assume that 85% of all SGRBs are created by the merger of two compact objects, the factor f_SGR = 0.85 yields a conservative limit on the merger rate. In order to create a relativistic outflow, a torus must be created around the newly formed BH. Simulations and theoretical analyses show that the formation of a GRB depends on several parameters, such as the spin of the BH, the mass ratio of the binary, or the compactness of the NS (Pannarale et al. 2011; Rezzolla et al. 2011). Since these parameters are hard to generalize, we simply assume that all compact object mergers produce a GRB. GRBs are believed to emit their radiation in a collimated cone. The half-opening angle θ defines the fraction of the sky where the burst can be seen, f_b = 1 − cos θ. The angle θ is highly uncertain, especially in the case of SGRBs, as it depends on the model and the Lorentz factor of the outflow (Panaitescu & Kumar 2001). Measurements of SGRB half-opening angles range from a few degrees to over 25° (Soderberg et al. 2004; Burrows et al. 2006; Grupe et al. 2006; Panaitescu 2006; Racusin et al. 2009). In the following, we set 1/f_b = 15, corresponding to a half-opening angle θ ≃ 20° (Bartos et al. 2011). By taking into account all of these factors, the merger rate R_merger(z) and the expected rate of SGRB observations R_SGRB(z) are

$$\begin{equation} R_{\mathrm{merger}}(z)=\frac{f_\mathrm{SGR}}{T\,f_b\,f_{{\rm FOV}}\,f_R} N^{\prime }(z) \end{equation} \tag{ 6 }$$

and

$$\begin{equation} R_{\mathrm{SGRB}}(z)=\frac{1}{T} N^{\prime }(z)\,, \end{equation} \tag{ 7 }$$

respectively.

The results of our analysis are summarized in Table 2. The approximate local merger rate of binary compact objects ranges from 479 to 1025 Gpc⁻³ yr⁻¹, depending on the chosen rate function. The upper limit comes from the model with unity rate function, which assumes no evolution on star formation over cosmological distances. The Porciani delayed rate function with delay times of 20 Myr and 100 Myr gives merger rates about 50% larger than without the delay. Extrapolating this result to the other functions, a reasonable estimate for the merger rate in the local universe is in the range ≃ 500–1500 Gpc⁻³ yr⁻¹. These results strongly depend on the half-opening angle of the SGRB jets. Figure 3 shows the dependence of the merger rate on the opening angle θ. The smaller the angle, the larger the number of mergers because the observer must be in the outflow cone to detect the SGRB. Assuming a half-opening angle of 10°, the merger rate could be as high as several thousand Mpc⁻³ yr⁻¹. A more isotropic large half-opening angle of 60° yields a rate of the order of 100 Mpc⁻³ yr⁻¹.

Zoom In Zoom Out Reset image size

Assuming that a satellite with a field of view comparable to the field of view of SWIFT is operating at the time of advanced GW detectors, and using the expected range for Advanced LIGO/Virgo (LIGO Scientific Collaboration and Virgo Collaboration 2010), we can estimate the number of coincident observations of SGRBs with GW counterparts. Under the above assumptions, we estimate about 0.2–1 coincident observations per year for an NS–NS merger progenitor (detector range ≃ 450 Mpc) and about 1–3 coincident observations per year for an NS–BH progenitor (detector range ≃ 930 Mpc), a result consistent with earlier estimates (see Metzger & Berger 2012; Bartos et al. 2012 and references therein). These values include the systematic uncertainties underlying the estimates, as explained above. Since the advanced detectors are expected to operate for several years, a few observations of coincident SGRB-GW events seem likely.

These estimates could improve significantly with a network of operating GRB satellites. Assuming that Fermi⁴ will be operating during the advanced detector era, as well as the planned SVOM mission (Paul et al. 2011) and Lobster,⁵ the coincident SGRB-GW detection rate could be higher than the above estimate by a factor ≃ 3. Ensuring that at least one GRB mission is operational at the time of advanced detectors will be crucial for identifying the host galaxy, measuring its redshift and star formation rate, and gaining valuable astrophysical information.

Our estimates can be compared to earlier published results that were obtained with alternative approaches. The two most common methods that are used to estimate the merger rate of compact objects rely on deriving the rates from observed pulsar observations (Kalogera et al. 2004) or employing population synthesis models (Belczynski et al. 2007). Both approaches inherit large statistical or systematical errors. A recent review article summarizes these results, concluding that the rate of merger events is somewhere between 10 and 10000 Gpc⁻³ yr⁻¹ (LIGO Scientific Collaboration and Virgo Collaboration 2010). Other investigations rely on methods that are more similar to the method used here. Guetta & Piran (2005) use a sample of five SGRBs to estimate a merger rate in the range 8–30 Gpc⁻³ yr⁻¹. An earlier analysis by one of the authors, which is based on a less restrictive data sample and neglects individual GRB luminosities, yields a much higher rate of about 7800 Gpc⁻³ yr⁻¹ (Dietz 2011). Finally, a recent study by Coward et al. (2012) uses a different method based on single GRB observations. In this approach, the maximum distance at which individual SGRBs can be detected by SWIFT is calculated and the results are then combined to estimate a final local rate of 0.16–1100 Gpc⁻³ yr⁻¹. These results show that our estimates are consistent with, and confirm, previous merger rate estimates. Future SGRB observations and improved statistics may further strengthen this conclusion.

This work is the result of a Research Experience for Undergraduates (REU) project by Carlo Enrico Petrillo at the University of Mississippi. C.E.P., A.D., and M.C. are partially supported by the National Science Foundation through awards PHY-0757937 and PHY-1067985. The authors thank Jocelyn Read, Emanuele Berti, Maurizio Paolillo, Neil Gehrels, and Richard O'Shaughnessy for their help and valuable comments. This publication has been assigned LIGO Document Number LIGO-P1200015.

This Appendix describes the various rate functions that are used in the analysis.

Porciani. The Porciani rate function is the SF2 function in Porciani & Madau (2001):

$$\begin{equation} R(z) \propto \frac{\exp {(3.4\,z)}}{\exp {(3.4\,z)}+22}\,. \end{equation} \tag{ A1 }$$

Hernquist. The Hernquist rate function is (Hernquist & Springel 2003)

$$\begin{equation} R(z) \propto \frac{\chi ^2}{1+\alpha (\chi -1)^3\exp {(\beta \chi ^{7/4})}}\,, \end{equation} \tag{ A2 }$$

where

$$\begin{equation} H(z) = H_0\,\sqrt{(1+z)^3\Omega _M + \Omega _\Lambda }\,, \chi (z) = \left[\frac{H(z)}{H_0}\right]^{2/3}\,, \end{equation} \tag{ A3 }$$

and α = 0.012, β = 0.041.

Fardal. The Fardal rate function is defined as (Fardal et al. 2007)

$$\begin{equation} R(z) \propto \frac{a^{-p_2}}{(1+p_1\,a^{-p_2})^{p_3+1}} H(z), \end{equation} \tag{ A4 }$$

where p₁ = 0.075, p₂ = 3.7, p₃ = 0.84, and a = (1 + z)⁻¹.

Cole. The Cole rate function is

$$\begin{equation} R(z) \propto \frac{a+bz}{1+(z/c)^d}H(z), \end{equation} \tag{ A5 }$$

where the parameters are summarized for different authors in Table 3.

Delayed functions. Since the time between the formation of the compact objects and the merger of the binary system is typically of the order of the Gyr, the merger rate may not follow directly from the star formation rate. Assuming a delay with respect to the star formation rate, the delayed rate function is defined as

$$\begin{equation} R_t(z) = \int _0^{t_d} dt \frac{1}{1+z_f} R(z_{{\rm ret}})\,P(t)\,, \end{equation} \tag{ A6 }$$

where P(t) represents the probability distribution of the delay time. Population synthesis models (Belczynski et al. 2001; Postnov & Yungelson 2005) suggest that this distribution is a power law

$$\begin{equation} P(t)\propto t^\alpha \,, \end{equation} \tag{ A7 }$$

where α ≃ −1 for t > t_tmin. Although the observational literature has applied a much broader range of functional forms, a time delay probability distribution P(t) ∼ 1/t is sufficient for the purposes of this analysis, where binaries are produced in the field.⁶ The retarded redshift, i.e., the redshift at the time when the compact objects are formed, is

$$\begin{equation} z_\mathrm{ret}=T^{-1} \left(T(z)+t \right). \end{equation} \tag{ A8 }$$

The look-back time at redshift z is (see, e.g., Hogg 1999)

$$\begin{equation} T(z) = T_0 \int _0^{dz^{\prime }} \frac{dz^{\prime }}{(1+z^{\prime })\sqrt{(\Omega _k(1+z^{\prime })^3 + \Omega _k(1+z^{\prime })^2 + \Omega _\Lambda)}}\,. \end{equation} \tag{ A9 }$$

The expression in Equation (A8) requires an integration and function inversion that generally need to be evaluated numerically. However, if flat cosmological models with Ω_k = 0 are considered, then it is possible to obtain an analytic expression for the look-back time. Substituting

$$\begin{equation} v = \Omega _k(1+z)^3 \end{equation} \tag{ A10 }$$

into Equation (A8), the integral takes the form

$$\begin{equation} T(z) = \frac{T_0}{3} \int _{v(0)}^{v(z)} dv^{\prime } \frac{1}{v\sqrt{v+\Omega _\Lambda }}\,. \end{equation} \tag{ A11 }$$

Integrating, it follows

$$\begin{equation} T(z) = \frac{T_0}{3\sqrt{\Omega _\Lambda }} \left(L(v_0)-L(v_1)\right)\,, \end{equation} \tag{ A12 }$$

where

$$\begin{equation} L(v) = \ln \left(\frac{\sqrt{v+\Omega _\Lambda }+\sqrt{\Omega _\Lambda }}{\sqrt{v+\Omega _\Lambda }-\sqrt{\Omega _\Lambda }} \right) \,. \end{equation} \tag{ A13 }$$

This analytic expression can be used to calculate the look-back time for a given redshift. Solving Equation (A12) for L(v₁), we find

$$\begin{equation} L(v_1)=L(v_0)-3\sqrt{\Omega _\Lambda }\frac{T}{T_H} \equiv \ln {E(T)}\,. \end{equation} \tag{ A14 }$$

The value of L(v₀) does not depend on z or T. Equation (A14) is a function of Ω_M and Ω_Λ, as one can see from Equations (A13) and (A10). Solving Equation (A13) for z, one finally obtains

$$\begin{equation} z(T) = \left(\frac{\Omega _\Lambda }{\Omega _M}\right)^{1/3} \left[ \left(\frac{1+E(T)}{1-E(T)}\right)^2 -1\right]^{1/3}-1\,. \end{equation} \tag{ A15 }$$