Cosmic Evolution of Long Gamma-Ray Burst Luminosity Can-Min Deng1,2, Xiang-GaoWang 1,2, Bei-Bei Guo1,2, Rui-JingLu1,2, Yuan-ZhuWang 1,2, Jun-JieWei 3,4,Xue-Feng Wu3,4, andEn-Wei Liang 1,2 6 1 0 2 Received ; accepted n a J 8 2 TheAstrophysicalJournal,000:000-000201500 00 ] E H . h p - o r t s a [ 1 v 5 4 6 7 0 . 1 0 6 1 : v i X r a 1GXU-NAOC Center for Astrophysics and Space Sciences, Department of Physics, Guangxi University,Nanning530004,China;[email protected];[email protected] 2GuangxiKeyLaboratory fortheRelativisticAstrophysics,Nanning530004,China 3PurpleMountainObservatory,ChineseAcademy ofSciences, Nanjing210008,China 4GraduateUniversityofChineseAcademyofSciences, Beijing100049,China – 2 – ABSTRACT Thecosmicevolutionofgamma-rayburst(GRB)luminosityisessentialforreveal- ing the GRB physics and for using GRBs as cosmological probes. We investigatethe luminosityevolution of long GRBs with a large sample of 258 Swift/BAT GRBs. Pa- rameterizedthepeakluminosityofindividualGRBsevolvesasL ∝ (1+z)k,weget p k = 1.49±0.19usingthenon-parametricτ statisticsmethodwithoutconsideringob- servationalbiasesofGRBtriggerandredshiftmeasurement. Bymodelingthesebiases with the observed peak flux and characterizing the peak luminosity function of long GRBs as a smoothly broken power-law with a break that evolves as L ∝ (1+ z)kb, b weobtaink = 1.14+0.99 throughsimulationsbasedonassumptionthatthelongGRB b −0.47 rate follows the star formation rate (SFR) incorporating with cosmic metallicity his- tory. The derived k and k values are systematicallysmallerthan that reported in pre- b vious papers. By removing the observational biases of the GRB trigger and redshift measurement based on our simulation analysis, we generate mock complete samples of258and1000GRBstoexaminehowthesebiasesaffectsontheτ statisticsmethod. Wegetk = 0.94±0.14andk = 0.80±0.09forthetwosamples,indicatingthatthese observationalbiases maylead to overestimatethek value. Withthelargeuncertain of k derived from our simulation analysis, one even cannot convincinglyargue a robust b evolutionfeatureoftheGRBluminosity. Subjectheadings: Gamma-rayburst: general–methods: statistical–stars:formation – 3 – 1. Introduction Gamma-ray bursts (GRBs) are the most luminous events in the universe. They have been detected up to a redshift of 9.4 so far (GRB 090429B; Cucchiara et al. 2011). It is generally believed that these events are from mergers of compact stellar binaries (Type I) or core collapses of massivestars (Type II) (e.g. Woosley 1993; Paczynski 1998; Woosley & Bloom 2006; Zhang 2006, Zhang et al. et al. 2007, 2009). The durations of prompt gamma-rays of Type I GRBs are usuallyshortandtheyarelongforTypeIIGRBs. Thismayleadtoatentativebimodaldistribution of the burst duration and two groups of GRBs were proposed, i.e., short GRBs (SGRBs) with T < 2 s and long GRBs (LGRBs) with T > 2 s, where T is the time interval of from 5% to 90 90 90 95% gamma-ray photons collected by a given instrument (Kouveliotou et al. 1993)1. Since the LGRB rate may follow the cosmic star formation rate (SFR), it has been proposed that LGRBs could be probes for the high-redshift universe, and may be used as a potential tracer of SFR at high redshift as where it becomes difficult for other methods (Totani 1997; Wijers et al. 1998; Blain & Natarajan 2000; Lamb & Reichart 2000; Porciani & Madau 2001; Piran 2004; Zhang & Me´sza´ros 2004; Zhang 2007; Robertson & Ellis 2012; Elliott et al. 2012; Wei et al. 2014; wang et al. 2015). It was also proposed that LGRBs may be promisingrulers to measure cosmological parameters and dark energy (e.g., Dai et al. 2004, Ghirlanda et al. 2004; Liang & Zhang et al. 2005). ThecosmicevolutionofGRBluminosityisessentialforrevealingtheGRBphysicsandusing GRBs as cosmological probes, but it is poorly known being due to lack of a large and complete GRB sample in redshift. With growing sample of high-redshift LGRBs since the launch of Swift 1Note that the statistical significance of the bimodal distribution depends on instruments and energybands(e.q.,Qinetal. 2013). Burstdurationonlyisdifficulttoclarifythephysicaloriginof some GRBs, such as GRBs 050724 and 060614 and new classification scheme and classification methodshavebeen proposed(Zhang 2006;Zhanget al. 2007,2009;Lu¨ et al. 2010). – 4 – satellite, it was found that GRBs rate increases significantly faster than the SFR, especially at high-z (Daigne et al. 2006; Le & Dermer 2007; Salvaterra & Chincarini 2007; Guetta & Piran 2007;Yu¨kseletal. 2008;Li2008;Salvaterraetal. 2008;Kistleretal. 2008). Itisunclearwhether this excess is due to some sort of evolutions in an intrinsic luminosity/massfunction (Salvaterra & Chincarini 2007; Salvaterra et al. 2009, 2012)or cosmicevolutionof GRB rate (e.g., Butler et al. 2010; Qin et al. 2010; Wanderman & Prian 2010; Wang & Dai 2011). Coward et al. (2013) even proposed that it is not necessary to invoke luminosityevolution with redshift to explain the observedGRB rate at high-zby carefully takingselection effects intoaccount (see also Howell& Coward, 2013) Attempts to measure the cosmic evolution of LGRB luminosityhave been made by several groups. The non-parametric τ statistical method (Efron & Petrosian 1992) is usually used to estimate the possible intrinsic L − z correlation by simplifying the evolution feature as L = L (1+z)k for individualGRBs, where the L is the luminosityof local GRBs. Strong p p,0 p,0 L − z dependence with a k value varying from 1.4 to 2.7 with different samples is reported (Lloyd-Ronning et al. 2002; Yonetoku et al. 2004; Kocevski & Liang 2006; Petrosian et al. 2009). Alternately, the comic evolution of GRB luminosity was also estimated with the cosmic evolution of their luminosity function. With this approach, the luminosity function is usually adopted as a broken power-law with a break at L , which evolves as L = L (1+z)kb. k then b b b,0 b can beestimatedby fittingtheobservedL and z distributionsviaMonteCarlo simulationsunder p assumption that GRB rate follows the SFR. Using this approach, the k derived from different b samplesisaround 2 (Salvaterraet al. 2012;Tan et al. 2013). Sample selection and observational biases are critical for measuring the GRB luminosity evolution feature. As mentioned above, results derived from different samples are significantly different. The redshift-known GRB samples are inevitably suffered observational biases on the flux truncation, trigger probability, redshift measurement, etc (e.g., Coward et al. 2013). By – 5 – modeling these biases with a large and uniform sample of GRBs with redshift measurement in a broad redshift range is essential to robustly estimate the luminosity evolution. The Swift missionhas establishedaconsiderablesampleof LGRBs withredshift measurementin arange of 0.1 ∼ 9.4 over 10 operation years. This paper revisits the cosmic luminosity evolution with the current redshift-known GRB sample by carefully considering the possible observational biases. We report our sample and the apparent luminosity dependence to redshift in §2. We derive the L−z dependence for GRBs from the current sample of Swift GRBs with redshift measurement with the non-parametric τ-statistics method (§3). By considering the BAT trigger probability and redshift measurement probability, we evaluate the L − z dependence with simulations by assuming that GRB rate follows the SFR incorporating with the cosmic metallicity history (§4). Based on oursimulationanalysis,we furtherinvestigatehowtheseobservationalbiases affect the results of the τ statistics method (§5). We present our conclusions and make brief discussion in §6. Throughout the paper the cosmological parameters H = 71 km s−1 Mpc−1, Ω = 0.3 and 0 M Ω = 0.7are adopted. Λ 2. Sampleselectionand Data Oursampleincludesonlytheredshift-knownlongGRBs observedwithSwift/BATfromJan. 2005toApril2015. Low-luminosityGRBsareexcludedsincetheymaybelongtoanotherdistinct population(e.g., Liang et al. 2007; Chapman et al. 2007; Lu¨ et al. 2010). Althoughthe durations of some GRBs, such as GRBs 050724 and 060614, are larger than 2 seconds, we do not include them in our samplesince they may be from merger of compact stars (Berger et al. 2005; Gehrels et al. 2006; Zhang et al. 2007; 2009). The afterglow and host galaxy observations of the high-z short GRB 090426(T = 1.2s) showthat it may be from collapseof a massivestar(Antonelli et 90 al. 2009;Xinet al. 2010). WethusincludethisGRB in oursample. Wefinally obtainasampleof 258 GRBs. They are reported in Table 1. Their redshifts (z), peak photonfluxes (P), and photon – 6 – indices(Γ)are taken fromtheNASAwebsite2. GRB radiation spectra are very broad. They can be well fitted with the so-called Band function, which is a smooth broken power-law characterized with low and high photonindices at a break energy (Band et al. 1993). Thepeak energy (E )of theνf spectrumranges from several p ν keVs to MeVs (Preece et al. 2000; Liang & Dai 2004; LU¨ et al. 2010; Zhang et al. 2011; von Kienlin et al. 2014). BAT energy band covers a very narrow range, i.e., 15-150 keV. The GRB spectraobservedwithBAT are adequatelyfitted withasinglepower-law,N(E) ∝ E−Γ (Zhang et al. 2007;Sakamotoet al. 2009). Sincelack ofE information,wedonot calculatethebolometric p luminosity of the GRBs in order to avoid unreasonable extrapolation with spectral information derivedfromBATdata. Wecalculatetheirpeakluminosity(L )usingthe1speakfluxesmeasured p in the15-150 keV band for ouranalysis. Thedistributionof theGRBs in thelogL −log(1+z) p plane is shown in Figure 1. A tentativecorrelation with large scatter is observed. The Spearman correlation analysis yieldsa correlation efficient of0.72 and a chance probabilityp < 10−4. With the ordinary least squares bisector algorithmand consideran intrinsicscatter of σ (e.g., Weiner in et al. 2006), we derive a relation of logL = (49.98± 0.09)+ (2.95±0.19)log[(1 + z)] and p σ = 0.72±0.02by takingtheluminosityuncertainty intoaccount3. in The flux thresholdof Swift/BAT is complicated, and the triggerprobability ofa GRB with a peak flux closeto theinstrumentthresholdis muchlowerthan that ofhigh-flux GRBs (e.g., Stern etal. 2001Qinetal. 2010). Inaddition,intheimagetriggermode,aGRBtriggeralsodependson theburst duration, hence theburst fluence (Band 2006, Sakamoto et al. 2007, Virgili et al. 2009). We adopt a flux threshold as 1.0×10−8 erg cm−2 s−1, which is availableby the BAT team in the 2http://swift.gsfc.nasa.gov/ 3Since the errors are very small and even no uncertainty is reported for the redshifts of most GRBs, wedo nottakeintoaccount theuncertaintiesofredshiftsin ourfit. – 7 – NASA webpage4.It isalsoshowninFigure1. 3. Measuring the intrinsicL −z dependence withthe τ statistics method p Because thesampleis greatly suffered from flux truncationeffect, it is uncertain whetherthe apparent L −z relation is intrinsic or only due to the truncation effect. The non-parametric τ p statisticalmethodisasimple,straightforwardwaytoestimatethetruncationeffectonanobserved correlation (Efron & Petrosian 1992; see also Lloyd-Ronning et al. 2002, Yonetoku et al. 2004 and Kocevski & Liang 2006). We first apply this method to measure the intrinsic dependence of L to(1+z) foroursample. We outlinethistechniqueas following. p Firstly,wepickup an associatedset of(L , z ), i i J = {j |L > L , z < z }, for 1 ≤ i ≤ 258 , (1) i j i j i,lim where z is the redshift corresponding to the flux threshold for a GRB of (L ,z ), and N is the i,lim i i size ofour sample. Thenumberof theset {J } is defined as N . Independence between L and z i i i i wouldmakethenumberofthesample R = number{j ∈ J |z ≤ z }, (2) i i j i uniformly distribute between 1 and N . To quantify the correlation degree between the two i quantities,theteststatisticτ isintroducedas (R −E ) τ = i i i , (3) P V i i pP whereE = (N +1)/2andV = (N2−1)/12aretheexpectedmeanandvariancefortheuniform i i i i distribution,respectively. Thisτ valueisemployedtomeasurethedatacorrelationdegreebetween 4 http://swift.gsfc.nasa.gov/about swift/bat desc.html – 8 – the two quantities of L and z in units of standard deviation. In case of that the two quantities are intrinsicallyindependentwithoutanycorrelation,theτ statisticsgivesτ = 0. Notethat N should i be greater than 5 since this method assumes Gaussian statistics. For the details of this technique pleaserefer toEfron &Petrosian (1992). We simply parameterize the cosmic evolution of the GRB luminosity as L = L (1 + z)k 0 for individual GRBs (see also Lloyd-Ronning et al. 2002; Yonetoku et al. 2004). We vary the k value and derive the corresponding τ value. Figure 2 shows τ as a function of k for our GRB sample. It is found that |τ| = 7.7 at k = 0, indicating that the null hypothesis that L and z are statisticallyindependent is rejected in a confidence levelof 7.7σ. In the case of k = 1.49±0.19, we have τ = 0 within 1 σ significance. Therefore, the intrinsic dependence of L to 1 + z is p L ∝ (1+z)1.49±0.19. It is muchshallowerthan theapparent one, which is L ∝ (1+z)2.95±0.19. p p Theapparent logL −log(1+z) relation isdominatedby theinstrumentselectioneffect. p 4. Measuring the intrinsicL −z dependence withthe comicsevolutionoftheluminosity p function by Monte Carlosimulations The evolutionof GRB luminositycan also showsup as the evolutionof theGRB luminosity function. The GRB luminosity function is usually described with a smooth broken power-law (e.g., Lianget al. 2007and references therein), L α1 L α2 −1 Φ(L ) = Φ + , (4) p p,0 (cid:20)(cid:18)L (cid:19) (cid:18)L (cid:19) (cid:21) b b where Φ is a normalization constant, α and α are the low and high indices of the luminosity p,0 1 2 function. We suggestthat theevolutionofGRB luminosityis parameterized with theevolutionof – 9 – L intheformof5 b L = L (1+z)kb, (5) b b,0 whereL istheL at z = 0. b,0 b We first show how L may evolve with redshift from the observational data. We plot the b cumulativeluminositydistributionsin different redshift ranges in Figure 3. One can observethat they show as broken power-laws, with a larger break luminosity (L ) at higher redshift range. b We fit them with Eq. 4. The results are reported in Table 2. One can find that L shows a trend b that GRBs at higher redshift tends to be more luminous, but one still cannot statistically claim a correlation between the break luminosity and the redshift based on the current sample since the uncertainty of the L is very large, especially in the redshift intervals of above z > 1. In the b redshift intervalof0 < z < 1, L is smallerthan thatin thehigh redshift intervalswithabout one b orderofmagnitude. However,weshouldnotethat theGRB subsetin thisredshiftintervalmay be contaminatedbythelowluminosityGRB populationas wementionabove. We measure the k valueby fitting the observed L and z distributionsin a physical context b p that the GRB rate follows the SFR rate incorporating the cosmic metallicityhistory by modeling theobservationalbiasesas reported inQin et al. (2010)viaMonteCarlo simulations. We assumethat theGRB rate(R ) traces theSFR and metallicityhistory(e.g., Kistleret LGRB al. 2008; Li 2008; Qin et al. 2010; Virgiliet al. 2011; Lu et al. 2012). TheintrinsicGRB redshift distributionthen can bewrittenas SFR(z)Θ(ǫ,z)dV(z) ρ(z) ∝ , (6) 1+z dz 5Physically, the luminosity evolution parameterized as L = L (1 + z)k for individual GRBs 0 is comparable to that parameterized as L = L (1 + z)kb for the break of the GRB luminosity b b,0 function,butthederivedk and k valuesfrom agivensamplewouldnotbecompletelythesame. b – 10 – where R (z) is the star formation rate (taken the form parameterized by Yu¨ksel et al 2008), SFR Θ(ǫ,z)isthefractionalmassdensitybelongingtometallicitybelowǫZ atagivenredshiftz(Z J J is the solar metal abundance) and ǫ is determined by the metallicity threshold for the production of LGRBs (e.g. Langer & Norman 2006). dV(z) is the co-movingvolumeelement at redshift z in dz aflat Λ CDM universe,givenby dV(z) c 4πD2(z) = L , (7) dz H (1+z)2[Ω (1+z)3 +Ω ]1/2 0 M Λ whereD (z) istheluminositydistanceatz and c isthespeedoflight. L We make Monte Carlo simulations based on Eqs. (4)–(7). The parameters of our model include α ,α ,L ,k and ǫ. We set these parameters in broad ranges, i.e., α ∈ (0,2), 1 2 b,0 b 1 α ∈ (1.8,4.5), L ∈ (1050,1053) erg/s, k ∈ (−0.5,3.0) and ǫ ∈ (0,1.6). The details of our 2 b,0 b simulation method please refer to Qin et al. (2010). We outline the simulation procedure as following. • We randomly pick up an ǫ value and generate a redshift z from the probability distribution ofEq. (6),thenrandomlypickupasetofluminosityfunctionparameters{L ,k ,α ,α } b,0 b 1 2 tocalculate theprobabilitydistributionsofL at z withEqs. (4)and (5). We then randomly p pickup aL atredhsift z. A simulatedGRB is characterized with asetof(L , z). p p • We calculatethepeak flux (F )in 15-150 keV band fora simulatedGRB, assumingaGRB p photon index as Γ = 1.69, which is the average value of the GRB in our BAT sample. The peak flux is compared withthe BAT flux thresholdto evaluatewhetherit can be observable withBAT. • Wesimulatethetriggerprobabilityandredshiftmeasurementprobabilitybymodelingthese probabilitiesas a function of the peak flux. The details please see Eqs. (9) and (10) of Qin et al. (2010).