Merger delay time distribution of extended emission short GRBs

The most popular progenitor model for short duration Gamma-Ray bursts (sGRBs) is the merger of two compact objects. However, the short GRB population exhibit a certain diversity: some bursts display an extended emission (EE), continuing in soft $\gamma$-rays for a few hundreds of seconds post the initial short pulse. It is currently unclear whether the origin of such bursts is linked to compact object mergers. Within the merger hypothesis, the redshift ($z$) distribution of short GRBs is influenced by the merger delay time, i.e., time elapsed between the merger and the formation of the binary star system, which is dominated by the time-scale for gravitational wave losses during the compact binary phase. We examine redshift distributions of short GRBs with extended emission to see whether their formation channel requires considerable delay post the star formation episode. Our results show that the $z$ distribution of EE bursts is consistent with the merger model. We attempted to compare the delay time distribution of the EE and the non-EE short bursts. However, no statistically significant difference could be seen within the limited sample size.


INTRODUCTION
The bimodal distribution of Gamma Ray Bursts (GRBs) in their temporal-spectral plane, as long-soft and short-hard bursts (Kouveliotou et al. 1993), calls for a distinction in their progenitor models. The most popular model of short duration Gamma Ray Bursts is the merger of compact binaries under loss of orbital energy and angular momentum through Gravitational Waves (GWs). Leading merger models of short bursts involve double neutron star (DNS) and Neutron Star Black-Hole (NSBH) systems (Paczynski 1986;Eichler et al. 1989).
A direct confirmation of the merger hypothesis has to wait till a joint detection by gravitational wave and electromagnetic detectors. However, several indirect observational characteristics of short bursts are favourable to this hypothesis. Long GRBs show a conclusive association with corecollapse supernovae (Galama et al. 1998;Hjorth et al. 2003), while short bursts do not. Positional offsets from the photo center of host galaxies are systematically larger than that of long bursts or core collapse supernovae, indicating natal kick velocities typical of neutron stars that can propel the binary off its birth place (Fong & Berger 2013). sGRB afterglows, being systematically fainter than that of long ones (Fong et al. 2015), may demand a lower density ambient medium typical of high galactic latitudes where compact binaries are expected to reach due to natal kicks. Short GRB hosts are of diverse types, implying a broad distribution in the time elapsed since the star formation epoch . Finally, recent claims of infra-red excess in short GRB afterglows is a promising signature of merger novae, arising from tidally disrupted neutron rich material from the merger (Li & Chevalier 1999;Tanvir et al. 2013).
Under the merger hypothesis, the redshift (z) distribution of sGRBs is a convolution of the star formation history of the universe and the merger delay time, where the latter is the time elapsed between the formation of a stellar binary and its eventual merger as two compact objects. The delay time is divided into two phases. First, the time taken for the two stars to evolve and become compact objects, which is typically of the order of 10 6 yrs. A major fraction of the delay time, ∼ 10 9 yrs (Hulse & Taylor 1975;Postnov & Yungelson 2014), is spent in the second phase, as two compact objects spiralling in by the emission of gravitational waves. Hence, merger delay time is an indication of the gravitational wave loss time-scale, and the z-distribution will be carrying an imprint of this.
As mentioned earlier, there are two progenitor channels, DNS and NS-BH systems, typically discussed in the context of short GRBs. Even within the DNS scenario, the central engine of the burst, formed as the result of the merger, can be diverse: a BH-torus system, a hypermassive short lived NS, or a long-lived NS, depending on the mass and equation of state of the components. Hence, it is likely that there could be a diversity within the short bursts themselves, and indeed some observations indicate the existence of such.
A unique characteristic shown by some short bursts is the presence of soft γ-ray emission (2 − 25 keV) extending to several hundreds of seconds post the initial short hard pulse (Villasenor et al. 2005). In some cases, the energy emitted in the extended emission is larger, even by a factor of 30 (Perley et al. 2009), than that in the prompt spike. Around 25% of short bursts are associated with extended emission (EE) (Norris et al. 2010). This long lasting prompt emission challenges models invoking a BH-torus system at the end of the merger. One of the alternate propositions is magnetar central engines, forming from DNS binaries or from accretion induced collapse of white dwarf binaries (Metzger et al. 2008).
The absence of spectral lag in the prompt emission, a characteristic differentiating short bursts from their long duration cousins, is shared by both EE and non-EE bursts (Norris & Bonnell 2006). However, there is no conclusive evidence for any defining characteristics the EE population may have from the remaining sGRB population. Though Troja et al. (2008) claimed that EE bursts potentially have a smaller host offset compared to other short bursts, studies based on a larger sample found that there are no conclusive differences in the offsets or host types of EE bursts from those of other short bursts (Fong & Berger 2013;Fong et al. 2013).
It is possible that EE-sGRBs are a distinct class which originates from a separate type of merging compact binaries compared to the rest of the short bursts. Or, the EE component could be an indication to a different kind of central engine. Since the timescale for merger under gravitational wave radiation losses depends on the component masses, eccentricity, and orbital period of the binary (see eq-3 in Postnov & Yungelson (2014)), naively one can expect different merger time distributions for different classes of merging binaries.
In this paper, we infer the delay time distribution of EE short GRBs from their redshift distribution under the ambit of the merger model. If EE short GRBs have a massive star origin similar to long GRBs, we expect to see its reflection in the inferred delay time distribution. Further, we compare the delay time of EE bursts with that of the non-EE short GRBs, albeit the caveat that the relative faintness of EE emission may lead to non-detections in higher redshifts. If a distinction in progenitor channel exists for the EE bursts, an imprint of the delay time would be visible in the redshift distribution in comparison with other short bursts. In section-2, we describe the sample selection. Theoretical modelling of the redshift distribution is explained in section-3. In section-4 we describe the results of our analysis. We summarize the paper in section-5.

SAMPLE OF SHORT GRBS
We first constructed a sample of short duration GRBs detected by Swift with reliable redshift measurements. To keep the sample uniform, we did not consider bursts detected by instruments other than Swift (for example, GRB 050709 detected by HETE is not part of the sample). This allows us to simplify the problem without having to consider varied biases across instruments. Since nearly all redshift measurements are available for Swift bursts alone, restricting to Swift does not compromise much on the sample size.
We used the 3rd Swift BAT catalogue (Lien et al. 2015) and the compilation by (Fong et al. 2015) to construct the sample. Both authors use T 90 < 2 sec criteria to categorize a burst as short. Our full GRB list contains only bursts up to 2015 October, where the BAT catalogue ends 1 . Redshift measurements are available only till March 2015 in Fong et al. (2015). For the one burst in our sample beyond this period, GRB150423A we used the BAT catalogue and the online table of Jochen Greiner 2 for z measurement. For GRB 080123, we use the z value given in (Leibler & Berger 2010), since (Fong et al. 2016) does not list this burst.
We excluded bursts with confusing redshift information (for example, 050813 and 051210) from the sample. GRBs 060505 and 090927, displaying spectral lag or hardness ratio typical of long bursts are not included. The SN-less 100 s duration GRB 060614 is also excluded. Two other bursts, 070506 and 090530, which are debated to be short bursts but having a main pulse longer than 2 s are also not part of the sample.

EE sample
(Lien et al. 2015) list short bursts with definite (table-3 of original paper) and possible (table-4 of original paper) extended emission. This sample (Lien-sample for the rest of the paper), is a collection of all previous literature including GCN circulars on EE searches. Possible shorts are the ones where the pulse is slightly longer than 2 s, or when the EE is weak, or if its existence is confusing due to stronger sources in the field of view.
We cross checked the list with the T 90 < 2 s bursts in five authors searching extended emission in GRBs (Gompertz et al. 2014;D'Avanzo et al. 2014;Kagawa et al. 2015;Kaneko et al. 2015;Fong et al. 2016). The sample selection criteria used by various authors are different, leading to a certain level of disagreement in results.
All T 90 < 2 s bursts for which Gompertz et al. (2014) detects EE are classified as EE in Lien et al. (2015). Their sample extends only till 2011. All T 90 < 2 s bursts of the DÁvanzo sample are included in the Lien-sample. Since they exclude bursts at low Galactic extinction lines of sight, their sample is much smaller than that of Lien. Kagawa finds three short bursts with EE which are also part of the Lien-sample. Eight of Kaneko's EE bursts have z information, of which the five with T 90 < 2 s are part of the Lien-sample. Duration of two bursts (051016B, 070506) are larger than 2 sec and the classification of 090927 is dubious. Hence, their absence in Lien-sample is justified.
GRB 061006 and 061210 are absent from the sample of Kagawa and DÁvanzo due to their selection criteria However, Lien et al. (2015) report EE in some short bursts which other authors failed to detect. Kagawa et al. (2015) . 2015) and (Fong et al. 2015). We have also considered several other literature on sGRB EE emission before arriving at the EE sample. See text for details. The last column indicates whether a burst is part of the EE sample or not.
Hence, we see that none of the conflicts between the authors are relevant here, and proceed to use Lien et al. (2015) as the reference for EE bursts. Our definite-EE sample has all Lien EE-bursts with confirmed redshift measurements (excluded the ones having only upper or lower limits of z), which is 5 bursts. For the full EE sample, we add two weak EE detections from Lien's table-4. Other bursts classified as possible EE in their table-4 are not considered because they have longer durations or the EE detection is confusing. In table-1, we lists all the bursts in our sample. Figure-1 is the cumulative z distribution of the bursts. It is obvious from the distribution that the EE bursts are confined to smaller redshifts. This could potentially come from a selection bias if EE bursts are fainter on an average, hampering their z determination. We have not included GRB160410A in this figure, which we include in a final special run of simulations alone. We proceed to model the cumulative z distribution of the EE bursts under the merger hypothesis to see the most agreeable merger delay time distribution. We also do a comparison between the EE and the non-EE sample. But the small size of the EE sample can affect the statistical significance of the result while comparing the EE and non-EE sample. Moreover a potential selection bias due to the faintness of the EE emission can interfere in the non-EE sample (if the EE is too faint to be detected, the burst will get classified as a non-EE burst) (Perley et al. 2009;Norris et al. 2010).

MODELLING REDSHIFT DISTRIBUTION
The broad framework we adopt to model the redshift distribution is well established in the literature (Guetta & Piran 2005;Nakar et al. 2006;Petrillo et al. 2013;Hao & Yuan 2013;Wanderman & Piran 2015). Nevertheless, for the sake of completion, we describe it below.
Rate of compact object mergers per unit co-moving volume, ρ merg , is a convolution of the star formation history of the universe ρ SFR (z) and delay time distribution f (τ) representing the probability that the merger happens with a time delay τ.
τ, the time elapsed between two epochs represented by z and z , can be written as (t(z)−t(z ) where t(z) is the age of the universe at redshift z. The coefficient A here takes care of the fraction of stars ending up as the kind of compact binaries in consideration. This fraction is assumed to be independent of redshift, and hence will not enter the normalized expression.
If all mergers lead to short GRBs, the rate N SHB of short hard bursts per unit z per unit time detected by an ideal detector is given by, However, for a detector with a limiting flux sensitivity of f lim which will not be able to detect bursts of all luminosity, the number of detectable bursts will depend on the luminosity function φ(L). Hence, where L min and L max are the lowest and highest luminosities in the burst distribution. If L min < L lim , the minimum detectable luminosity corresponding to f lim , L min will be replaced with L lim , which is a function of distance and hence z.
Hence, the normalized cumulative distribution at a redshift z 0 can be written as, Thus, the major components of the model are the star formation history (ρ SFR (z)), delay time distribution f (τ), and the luminosity function φ(L) of short GRBs. While there is a reasonable understanding of the cosmic star formation history upto a moderate redshift (Madau & Dickinson 2014), the merger delay time distribution and short GRB luminosity function are not well understood. For both these components we will have to resort to empirical prescriptions.

Components of the model
Star formation history of the universe: There is considerable uncertainity in the star formation history beyond z ∼ 2.
Madau & Dickinson (2014) uses a set of post-2006 galaxy surveys in rest-frame Far-UV, Mid, and Far IR to extract a ready-to-use functional form of the cosmic star formation history. The survey data they use is restricted to z < 8. The best fitting function is given by, ρ SFR (z) = 0.015 (1 + z) 2.7 1 + [(1 + z)/2.9] 5.6 M year −1 Mpc −3 .
This function rises in the early universe (between 3 z 8), peaks at a redshift of around 1.5 -2, and gradually declines to the present day universe. It shows a steady behaviour in the high-z range. There is some debate in the literature about the nature of the SFR for z > 9 (see Madau & Dickinson (2014) for details). Though such large redshifts are important to our analysis, since there is no definitive answer to the nature of high-z SFR, we resort to the (Madau & Dickinson 2014) prescription which uses the most recent state-of-the-art surveys. Delay time distributions: In the recipe for delay time, we are ignoring the time taken for the two main sequence stars to evolve into compact objects. The general expression for gravitational wave loss time scale for a compact binary depends on the period, eccentricity, the reduced mass, and the total mass of the binary system (Postnov & Yungelson 2014).
The simplest empirical form of delay time distribution f (τ), is a power-law function in τ. If the initial orbital separation a i of the binaries have a power-law distribution of the form f (a i ) ∝ a q i , the gravitational wave delay time will be distributed as f (τ) ∝ τ (q−3)/4 (Piran 1992). This calculation ignores any possible spread in orbital eccentricity. Assuming q = −1 in accordance with the X-ray binary data in our galaxy, Piran (1992)  Previous studies of short GRB rates in the literature have used simple empirical delay time distributions, mostly power-law (PL) and log-normal (LN) functions in τ (Nakar et al. 2006). The PL function is characterised by index η, while the LN function depends on two parameters, the mean delay time τ and the standard deviation σ. The LN model is an approximation of a constant delay time distribution. We also use these two empirical delay time functions in our calculations. We have used a lower cut-off of 20 Myr in the PL model, which represents the typical timescale for single star evolution. At a given z, we use the age of the universe for the corresponding z as an upper limit to τ. In figure-2 we show the sensitivity of the model N (< z 0 ) on η, τ , and σ. We use inferences from these dependences to constrain the parameter space (see next section).
Short GRB luminosity function: Luminosity function of short duration GRBs is convolved always with the rate, and manifests in the z distribution as well as in the peak flux distribution. Nakar et al. (2006) have used a single powerlaw luminosity function model, Φ(L) ∝ L −β , since the sample size was small during that era. They found the luminosity distribution is consistent with a β = 1.0. Wanderman & Piran (2015) found that a double power-law model is better fitted, with an index similar to (Nakar et al. 2006) at the lower end of the luminosity function. Since the sample of EE bursts are small, we chose to use the simplest form of the luminosity function, a single power-law throughout this work.
The minimum detectable luminosity at a distance d L (z) is, where k(z) is the cosmological k-correction which depends on the spectral shape of the burst. This factor appears in the equation since for a burst at redshift z, the BAT energy range (E 1 − E 2 ) corresponds to emitted photons of energy E 1 (1 + z) − E 2 (1 + z). Hence, the BAT flux of a burst with a cosmic rest-frame luminosity L will be proportional to (1 + z) 1−α where α is the spectral index ( f ν ∝ ν −α ). Most Swift bursts have an α ∼ 1 if fitted by single powerlaw spectral models (Lien et al. 2015). Hence we can safely ignore the (1+ z) factor in the equation-6. We use f lim = 5 × 10 −9 erg/s/cm 2 for the BAT threshold flux following Cao et al. (2011). For the range of z we encounter in this analysis, L lim varies from 10 45 to 10 51 and we always considered L min to be L lim . We restricted our analysis to β > 1 where L max remains irrelevant. Thus, eq-4 reduces to, Throughout this paper, we have used a standard cosmological model with H 0 = 69 km/s/Mpc, Ω m = 0.29, Ω λ = 0.71.

ANALYSIS AND RESULTS
Though the star formation rate drops down towards higher redshifts in our SFR model, we checked the convergence of the integration in eq-1. Under the assumption that the fitting function of Madau & Dickinson is rightly reproducing the unknown SFR beyond z of 8, we performed the integration for different values of z max and found that an z max ∼ 50 ensures convergence. We tested the model for (η > 10) indices of the PL delay time function and found that it is able to reduce the resulting ρ merg to the SFR.
For a given SFR, the parameter space reduces to the two dimensional (β, η) space for the PL model. We varied the power-law index η from −3. to 0.5. Since the variation in σ leads to a relatively lesser deviation between the models (see figure-2), we fixed σ at 0.2 Gyr and essentially reduced the parameters space again to the two dimensional (β, τ ) for the LN model. We increased σ to 0.6 Gyr but did not find any significant difference in N (< z 0 ). τ was made to vary from 2.5 to 4 Gyr, since a mean delay above 4 Gyr would have led to all bursts happening z < 1, which clearly is ruled out by the data. The index of the luminosity function β was varied from 1.1 to 3..
We constructed normalized cumulative distributions of both the data and the model. We normalized the model such that N (< z 0 ) = 1 at the highest redshift of the sample.
In order to infer the best fit delay time distribution, we undertook both χ 2 analysis and the Kolmogorov-Smrinov (KS) test. Since the sample is homogeneous, we did not go for a maximum likelihood estimate.
For the χ 2 analysis, we computed the summed square of the deviation between the data and the model. Both the tests were done on a grid of the parameter space.
There is a consistency between the KS and the χ 2 test results. The best fit values of the parameters from the χ 2 test is within the region of KS probability p > 0.05. Moreover, the region where χ 2 /dof ∼ 1 is well in agreement with the region of p > 0.05 in all the runs.

Delay time inferred for EE bursts
We find that the redshift distribution of the EE short GRBs is consistent within the predictions of the merger models (see figures 3-10). The EE bursts seem to originate after a considerable delay from the star formation episode. Under the log-normal delay time distribution model, a mean delay of ∼ 4 Gyr is required to produce the EE z distribution. Using a power-law model for the delay time distribution, the EE bursts seem to favour a flat to positive power-law index which indicates a major fraction with increased merger delay time (than a steep negative powerlaw). However, this could be due to the faintness of EE component which restricts their detection to larger redshifts. A distribution that stops at smaller redshifts could either be due to a longer delay time or due to faint emission. We included the highest z EE burst known to date, GRB160410A, to see if that can significantly reduce the delay time, but only for the PL model have we seen a change in model parameters. In PL model, as expected a trend towards a negative η is seen, which indicates that we may get more refined results with a larger sample of EE bursts.
As expected, the non-EE bursts in comparison are consistent with a lesser delay time. But again, the statistical significance of this result is not very high due to the potential selection bias. At high z values, EE need not be detected.
Our best fit luminosity function indices are well in agreement with that of the previous studies (Nakar et al. 2006). For both EE and non-EE samples, a β ∼ −2 is found to be consistent, which may be an indication of a single formation channel for both EE and non-EE bursts.
The major caveat of the analysis is the small sample size of EE bursts, especially them being restricted to relatively lower redshifts. The same can also contaminate the non-EE sample as some bursts might as well have a faint EE. Within this constraint, we infer that a delay time of the order of Gyrs is required in the formation of these bursts.
Since the sample of EE bursts with z information is small, we restricted our analysis to a preliminary form with simplified model functions. It is possible that the underlying luminosity function is more complex. A more robust analysis with more complicated model functions can be performed once the sample size of EE bursts increases.

SUMMARY
We model the redshift distribution of extended emission short GRBs under the ambit of the binary merger model to see the possible delay in the formation channel of this class of bursts. We find that the EE bursts require to have a considerable delay after the star formation episode in their formation channel indicating their origin in old stellar populations like compact stars. Luminosity function of both EE and non-EE bursts are consistent with each other, indicating their potential common origin. Current sample of short GRBs with z, especially of EE bursts are not sufficient to make a comparison of the delay time of bursts with and without extended emission. In future with a larger sample size, it will be possible to investigate their delay times separately.      , and could indicate that any difference in the merger delay time distribution of the EE sample w.r.t to the non-EE sample could merely be a selection effect. Figure 9. As a comparison, this figure shows the non-EE sample with the LN delay time model parameters. The best fit parameters shift to lower delay times compared to the EE sample, but the statistical significance of this difference is limited. Moreover, the somewhat higher delay time of the EE sample could also be due to the faintness that restricts it to lower redshifts. Figure 10. χ 2 values on grids and KS test contours for the PL delay time model for the non-EE sample. The colored region is where χ 2 ∼ DOF.