A MULTIVARIATE FIT LUMINOSITY FUNCTION AND WORLD MODEL FOR LONG GAMMA-RAY BURSTS

Ever since the discovery of the first gamma-ray burst (GRB) by the Vela satellites in 1967 (Klebesadel et al. 1973), there has been tremendous effort and attempts to constrain the energetics, luminosity function (LF), and the underlying mechanism responsible for these events. Early observations by the Konus (Mazets & Golenetskii 1981) and Ginga (Fenimore et al. 1988; Nishimura 1988) gamma/X-ray instruments suggested a possible link between GRBs and neutron stars with output energy ranges of the order of ∼10⁴⁰ erg. With the launch of the Compton Gamma-Ray Observatory (CGRO), the Burst And Transient Source Experiment (BATSE) on board CGRO dramatically changed the understanding of GRBs. While previous catalogs (e.g., Atteia et al. 1987) indicated an isotropic distribution of GRB sources, the BATSE observations extended this isotropy down to the weakest bursts. The non-homogenous (e.g., Fenimore et al. 1993) and isotropic spacial event distribution (e.g., Meegan et al. 1992; Briggs 1993; Fishman et al. 1994) provided, for the first time, strong support for a cosmological versus galactic origin of GRBs, undermining neutron stars in the local universe as the potential candidates for some—if not all—classes of gamma-ray events. Furthermore, the joint duration–hardness distribution of GRBs provided a direct evidence for at least two separate classes of GRBs: long-soft versus short-hard (e.g., Kouveliotou et al. 1993, see also Figure 1 here).

The possibility of a cosmological origin for GRBs indicated an enormous output energy on the order of ∼10⁵¹ erg (e.g., Dermer 1992). Nevertheless, an accurate description of the GRB LF also required knowledge of the GRB cosmic rate, information that could not be extracted from BATSE observations alone. This became possible only with the launch of the Italian-Dutch X-ray satellite BeppoSax (Boella et al. 1997) and the identification of the first GRB with firmly measured cosmological redshift (Metzger et al. 1997) that marked the beginning of the afterglow era in the field of GRBs. The launch of the Swift satellite (Gehrels et al. 2004) was another milestone that revolutionized the study of GRBs by facilitating the X-ray afterglow observations (Burrows et al. 2005) and further ground-based follow-ups for redshift measurement.

Alongside the observational triumphs over a few decades, several theoretical models have stood up against the rivals based on the available evidence and GRB data. Most prominently, the Collapsar model (e.g., Woosley 1993) has been relatively successful in linking the long-duration class of gamma-ray bursts (LGRBs) to the final stages in the lives of massive stars, while the short-duration class of bursts (SGRBs) is generally attributed to the coalescence of compact binary systems (e.g., Paczynski 1986; Nakar 2007 and references therein). The two classes of SGRBs and LGRBs in this work correspond to Type-I and Type-II GRBs, respectively, according to the physical classification scheme of Zhang et al. (2007) and Bloom et al. (2008). Further refinement of the potential candidates, as the progenitors and the emission mechanism for both classes, requires more rigorous analysis of observational data in all possible energy frequencies. In particular, the prompt gamma-ray emission of LGRBs has been subject of intense observational and theoretical studies.

Beginning with BATSE observations, numerous authors have examined the prompt emission of LGRBs searching for potential underlying correlations among the spectral parameters (e.g., Nemiroff et al. 1994; Fenimore et al. 1995; Mallozzi et al. 1995; Petrosian & Lee 1996; Brainerd 1997; Dezalay et al. 1997; Petrosian et al. 1999; Lloyd et al. 2000; Norris et al. 2005). The lack of known redshifts for BATSE events and poor knowledge of LGRB cosmic rates, however, strongly limited the prediction power of such analyses. Instead, the first direct evidence for potential correlations and constraints on the distributions of the prompt-emission spectral parameters came with a few LGRBs detected by BATSE, BeppoSax, IPN, or HETE-II satellites with measured redshifts (e.g., Reichart & Lamb 2001; Amati et al. 2002; Ghirlanda et al. 2004; Yonetoku et al. 2004) and was further developed by the inclusion of Swift LGRBs (e.g., Schaefer 2007; Gehrels et al. 2009). Such findings, however, have been criticized for relying primarily on a handful of events with high signal-to-noise ratio (S/N) required for spectral analysis with afterglows sufficiently bright for redshift measurement, arguing that the proposed joint distributions of the spectral properties do not represent the entire underlying population of LGRBs (e.g., Band & Preece 2005; Nakar & Piran 2005; Li 2007; Butler et al. 2007, 2009; Shahmoradi & Nemiroff 2009, 2010, 2011). Responding to criticisms, attempts were made to model the effects of different gamma-ray instruments' detection thresholds and the limiting effects of spectral analysis (e.g., Ghirlanda et al. 2008; Nava et al. 2008; see Shahmoradi & Nemiroff 2011 for a review of relevant literature).

Despite significant progress, difficulties in modeling the complex effects of detector threshold on the multivariate distribution of the prompt-emission spectral properties and the lack of a sufficiently large sample of uniformly detected LGRBs has led the community to focus on individual spectral variables, most importantly the LF (e.g., Petrosian 1993; Schmidt 1999, 2001, 2009; Kommers et al. 2000; Kumar & Piran 2000; Band 2001; Porciani & Madau 2001; Sethi & Bhargavi 2001; Stern et al. 2002; Guetta et al. 2005; Salvaterra & Chincarini 2007; Salvaterra et al. 2009; Campisi et al. 2010; Wanderman & Piran 2010; Salvaterra et al. 2012). Recently, Butler et al. (2010, hereafter B10) presented an elaborate multivariate analysis of Swift LGRBs, including the potential correlations among three temporal and spectral variables: total isotropic energy emission (E_iso (erg)), a definition of duration (T_r45 (s)), and the spectral peak energy (E_{p, z} (keV)) of the bursts. Focusing their analysis on the E_iso–E_{p, z}–T_r45 interrelation, B10 find a strong and significant correlation—but with a broad scatter—between the isotropic emission and the peak energy of the Swift LGRBs. Also realized by B10 is the possibility of a positive—and perhaps strong—correlation between the duration and the total energy output of the bursts. Moreover, to accommodate the potential existence of a large population of sub-luminous events, a broken power-law LF is used by B10.

It is known that the energy budget of GRBs must be limited and a turnover in the LF at the low-luminosity tail of the population is expected. However, the choice of the broken power law as the candidate LF is justified by the fact that current observational data cannot constrain this turnover point, which is expected to be far below the detection threshold of current gamma-ray instruments. The existence of a turnover point in the LF can have important clues for the underlying physics of LGRBs. The low-luminosity tail of the LF in such a case would be the result of the convolution of the stellar mass distribution with the LGRB rate as a function of the properties of massive dying stars, imposed by the mechanism.

Motivated by the search for the potential shape of the LGRB LF at low luminosity and the flurry of recent reports on the possible existence of correlations among LGRBs' temporal and spectral variables, the author of this paper has considered a wide variety of different statistical models that could incorporate and explain all observed correlations and spectral distributions, while preserving the prediction power of the model for the spectral properties of a potentially large fraction of LGRBs that could go undetected by current gamma-ray detectors. Here it will be shown that it is indeed possible to construct an LGRB world model capable of describing most (if not all) prompt-emission spectral and temporal properties observed in the current data catalogs. Toward this, the presented analysis is mainly focused on the largest catalog of GRBs available to date: the BATSE catalog of 2130 GRBs (Paciesas et al. 1999). The author shows that there is still an enormous amount of information buried in BATSE observations that need to be dug out by GRB researchers.

In the following sections, an example of such data mining on BATSE catalog will be presented: Section 2.1 presents an elaborate method for classifying GRBs into the two known subclasses of SGRBs and LGRBs. In Section 2.2 an LGRB world model capable of describing BATSE observations is presented, followed by a discussion of the procedure for fitting the model to data in Section 2.3. Univariate and multivariate goodness-of-fit (GoF) tests are performed to ensure that the model predicts BATSE data accurately (Section 2.4). The results from model fitting and implications on the multivariate distribution of the prompt-emission spectral properties are discussed in Section 3. A summary of the analysis and important outcomes of the proposed LGRB world model will be presented in Section 4.

2.1. Sample Selection

Depending on their detection criteria, some gamma-ray detectors might facilitate detection of one class of bursts over the others. For example, the specific detector sensitivity of the Burst Alert Telescope (BAT) on board the Swift satellite (Gehrels et al. 2004) results in better detections of LGRBs over SGRBs (e.g., Band 2003, 2006). Therefore, a simple GRB classification method, such as a cutoff in the observed duration distributions of GRBs (T₉₀ ∼ 3 s) generally results in LGRB samples that are minimally contaminated by SGRBs (e.g., B10). Compared to BAT, however, BATSE Large Area Detectors (LADs) had an increased relative sensitivity to SGRBs (e.g., Band 2003). Knowing that many temporal and spectral properties of LGRBs and SGRBs, most importantly T₉₀, overlap, identification of LGRBs in the BATSE catalog solely based on their T₉₀ will likely result in a significant number of misclassifications (e.g., Figure 1, center panel).

Zoom In Zoom Out Reset image size

Thus, to ensure correct analysis of a long-duration class of GRBs, it is first necessary to collect the least biased sample of events, all belonging to LGRB category. The word "bias" here refers to the bias that might be introduced when using the traditional definition of GRB classes, based on a sharp cutoff on the duration variable T₉₀ (Kouveliotou et al. 1993), as it is generally assumed by many GRB researchers (e.g., Guetta et al. 2005; B10; Campisi et al. 2010; Wanderman & Piran 2010) or not assumed or explicitly discussed (e.g., Salvaterra & Chincarini 2007; Salvaterra et al. 2009, 2012). This goal is achieved by testing extensive varieties of classification and clustering methods, most importantly, the fuzzy clustering algorithms of Rousseeuw et al. (1996) and the method of fuzzy C-means (Dunn 1973; Bezdek 1981). Each BATSE-catalog GRB is assigned a probability (i.e., class coefficient) of belonging to the LGRB (versus SGRB) population according to the choice of clustering algorithm and the set of GRB variables used. This can be any combination of the peak flux (P_50–300 photons s⁻¹ cm⁻²) in the BATSE detection energy range, in three different timescales: 64, 256, and 1024 ms; bolometric fluence (S_bol(erg cm⁻²)); the spectral peak energy (E_p (keV)) estimates by Shahmoradi & Nemiroff (2010); duration (T₉₀ (s)); and the fluence-to-peak-flux ratio (FPR (s)). Then GRBs with LGRB class coefficient >0.5 are flagged as long-duration class bursts. Overall, the fuzzy C-means classification method with the two GRB variables E_p and T₉₀ is preferred over other clustering methods and sets of GRB variables (cf. Appendix A). This leads to the selection of 1376 events as LGRBs out of the 1966 BATSE GRBs having measured temporal and spectral parameters mentioned above.¹

As a further safety check to ensure minimal contamination of the sample by SGRBs, the light curves of 291 bursts among 1966 BATSE GRBs with LGRB class coefficients in the range of 0.3–0.7 are visually inspected in the four main energy channels of BATSE LADs. This leads to reclassification of 17 events (originally flagged as LGRB by the clustering algorithm) to potentially SGRB or soft gamma repeater (SGR) events, and reclassification of seven events (originally flagged as SGRB by the clustering algorithm) to LGRBs. The result is a reduction in the size of the original LGRB sample from 1376 to 1366 (Table 1). It is notable that the inclusion of the uncertainties on the two GRB variables T₉₀ and E_p turns out to not have significant effects on the derived samples of the two GRB classes discussed above. Also, a classification based on T₅₀ instead of T₉₀ results in about the same samples for the two GRB classes with only a negligible difference of ∼0.7%.

Table 1. 1366 BATSE Catalog Triggers Classified as LGRBs

Trigger	Trigger	Trigger	Trigger	Trigger	Trigger	Trigger	Trigger	Trigger	Trigger	Trigger	Trigger	Trigger	Trigger
105	107	109	110	111	114	121	130	133	143	148	160	171	179
204	211	214	219	222	223	226	228	235	237	249	257	288	332
351	394	398	401	404	408	414	451	465	467	469	472	473	493
501	516	526	540	543	548	549	559	563	577	591	594	606	630
647	658	659	660	673	676	678	680	685	686	690	692	704	717
741	752	753	755	761	764	773	795	803	815	816	820	824	825
829	840	841	869	907	914	927	938	946	973	999	1009	1036	1039
1042	1046	1085	1086	1087	1114	1120	1122	1123	1125	1126	1141	1145	1148
1150	1152	1153	1156	1157	1159	1167	1190	1192	1196	1197	1200	1204	1213
1218	1221	1235	1244	1279	1288	1291	1298	1303	1306	1318	1382	1384	1385
1390	1396	1406	1416	1419	1425	1432	1439	1440	1446	1447	1449	1452	1456
1458	1467	1468	1472	1492	1515	1533	1540	1541	1551	1552	1558	1559	1561
1567	1574	1578	1579	1580	1586	1590	1601	1604	1606	1609	1611	1614	1623
1625	1626	1628	1642	1646	1651	1652	1653	1655	1656	1657	1660	1661	1663
1664	1667	1676	1687	1693	1700	1701	1704	1711	1712	1714	1717	1730	1731
1733	1734	1740	1742	1806	1807	1815	1819	1830	1883	1885	1886	1922	1924
1956	1967	1974	1982	1989	1993	1997	2018	2019	2035	2047	2053	2061	2067
2069	2070	2074	2077	2079	2080	2081	2083	2087	2090	2093	2101	2102	2105
2106	2110	2111	2112	2114	2119	2122	2123	2129	2133	2138	2140	2143	2148
2149	2151	2152	2156	2181	2187	2188	2189	2190	2191	2193	2197	2202	2203
2204	2207	2211	2213	2219	2228	2230	2232	2233	2240	2244	2252	2253	2254
2267	2276	2277	2287	2298	2304	2306	2309	2310	2311	2315	2316	2321	2324
2325	2328	2329	2340	2344	2345	2346	2347	2349	2362	2367	2371	2373	2375
2380	2381	2383	2385	2387	2391	2392	2393	2394	2405	2419	2423	2428	2429
2430	2432	2435	2436	2437	2438	2440	2441	2442	2443	2446	2447	2450	2451
2452	2453	2458	2460	2472	2476	2477	2482	2484	2495	2496	2500	2505	2508
2510	2511	2515	2519	2522	2528	2530	2533	2537	2541	2551	2560	2569	2570
2581	2586	2589	2593	2600	2603	2606	2608	2610	2611	2619	2620	2628	2634
2636	2640	2641	2660	2662	2663	2664	2665	2671	2677	2681	2688	2691	2695
2696	2697	2700	2703	2706	2709	2711	2719	2725	2727	2736	2749	2750	2751
2753	2767	2770	2774	2775	2780	2790	2793	2797	2798	2812	2815	2825	2830
2831	2843	2848	2850	2852	2853	2855	2856	2857	2862	2863	2864	2877	2880
2889	2890	2891	2897	2898	2900	2901	2913	2916	2917	2919	2922	2924	2925
2927	2929	2931	2932	2944	2945	2947	2948	2950	2951	2953	2958	2961	2980
2984	2985	2986	2990	2992	2993	2994	2996	2998	3001	3003	3005	3011	3012
3015	3017	3026	3028	3029	3032	3035	3040	3042	3055	3056	3057	3067	3068
3070	3071	3072	3074	3075	3076	3080	3084	3085	3088	3091	3093	3096	3100
3101	3102	3103	3105	3109	3110	3115	3119	3120	3127	3128	3129	3130	3131
3132	3134	3135	3136	3138	3139	3141	3142	3143	3153	3156	3159	3166	3167
3168	3171	3174	3177	3178	3193	3212	3217	3220	3227	3229	3237	3238	3241
3242	3245	3246	3247	3255	3256	3257	3259	3267	3269	3276	3279	3283	3284
3287	3290	3292	3301	3306	3307	3319	3320	3321	3322	3324	3330	3336	3339
3345	3347	3350	3351	3352	3356	3358	3364	3369	3370	3378	3403	3405	3407
3408	3415	3416	3436	3439	3448	3458	3465	3471	3472	3480	3481	3485	3486
3488	3489	3491	3493	3503	3505	3509	3511	3512	3514	3515	3516	3523	3527
3528	3552	3567	3569	3588	3593	3598	3608	3618	3634	3637	3648	3649	3654
3655	3658	3662	3663	3664	3671	3717	3733	3740	3745	3765	3766	3768	3771
3773	3776	3779	3788	3792	3800	3801	3805	3807	3811	3814	3815	3819	3840
3843	3853	3860	3869	3870	3871	3875	3879	3886	3890	3891	3892	3893	3899
3900	3901	3903	3905	3906	3908	3909	3912	3913	3914	3916	3917	3918	3924
3926	3929	3930	3935	3941	3954	4039	4048	4095	4146	4157	4216	4251	4312
4350	4368	4388	4556	4569	4653	4701	4710	4745	4814	4939	4959	5080	5255
5304	5305	5379	5387	5389	5407	5409	5411	5412	5415	5416	5417	5419	5420
5421	5423	5428	5429	5433	5434	5447	5450	5451	5454	5463	5464	5465	5466
5470	5472	5473	5474	5475	5476	5477	5478	5479	5480	5482	5483	5484	5486
5487	5489	5490	5492	5493	5494	5495	5497	5503	5504	5507	5508	5510	5512
5513	5515	5516	5517	5518	5523	5524	5526	5530	5531	5538	5539	5540	5541
5542	5545	5548	5551	5554	5555	5559	5563	5565	5566	5567	5569	5571	5572
5573	5574	5575	5581	5585	5589	5590	5591	5593	5594	5597	5601	5603	5604
5605	5606	5608	5610	5612	5614	5615	5617	5618	5621	5622	5624	5626	5627
5628	5632	5635	5637	5640	5644	5645	5646	5648	5654	5655	5667	5697	5704
5706	5713	5715	5716	5718	5719	5721	5723	5725	5726	5729	5731	5736	5773
5867	5890	5955	5983	5989	5995	6004	6082	6083	6090	6098	6100	6101	6102
6103	6104	6111	6113	6115	6118	6119	6124	6127	6128	6131	6137	6139	6141
6147	6151	6152	6154	6158	6159	6165	6167	6168	6176	6186	6188	6189	6190
6194	6198	6206	6222	6223	6225	6226	6227	6228	6233	6234	6241	6242	6243
6244	6249	6266	6267	6269	6270	6271	6272	6273	6274	6279	6280	6283	6285
6288	6295	6298	6300	6303	6304	6305	6306	6308	6309	6315	6317	6319	6320
6321	6322	6323	6328	6329	6330	6334	6335	6337	6339	6344	6345	6346	6349
6351	6353	6355	6369	6370	6375	6380	6388	6390	6395	6396	6397	6399	6400
6404	6405	6408	6409	6413	6414	6419	6422	6425	6435	6437	6440	6444	6446
6448	6450	6451	6453	6454	6472	6487	6489	6490	6498	6504	6519	6520	6521
6522	6523	6525	6528	6529	6531	6533	6534	6536	6538	6539	6544	6546	6550
6551	6552	6554	6557	6560	6564	6566	6576	6577	6578	6582	6583	6585	6587
6589	6590	6592	6593	6598	6600	6601	6602	6605	6610	6611	6613	6615	6616
6619	6620	6621	6622	6625	6629	6630	6631	6632	6642	6648	6649	6655	6657
6658	6665	6666	6670	6672	6673	6674	6676	6678	6683	6686	6694	6695	6698
6702	6707	6708	6720	6745	6762	6763	6764	6767	6774	6782	6796	6802	6814
6816	6830	6831	6853	6877	6880	6882	6884	6891	6892	6903	6911	6914	6917
6930	6935	6938	6963	6987	6989	7000	7012	7028	7030	7064	7087	7108	7110
7113	7116	7130	7147	7164	7167	7170	7172	7178	7183	7185	7191	7206	7207
7209	7213	7219	7228	7230	7247	7250	7255	7263	7285	7293	7295	7298	7301
7310	7318	7319	7322	7323	7328	7335	7343	7357	7358	7360	7369	7371	7374
7376	7377	7379	7381	7386	7387	7390	7403	7404	7429	7432	7433	7446	7451
7452	7457	7460	7464	7469	7475	7477	7481	7485	7486	7487	7488	7491	7493
7494	7497	7500	7502	7503	7504	7509	7515	7517	7518	7520	7523	7527	7528
7529	7532	7533	7535	7548	7549	7550	7551	7552	7560	7563	7564	7566	7567
7568	7573	7575	7576	7579	7580	7587	7588	7597	7598	7603	7604	7605	7606
7607	7608	7609	7614	7615	7617	7619	7625	7630	7635	7638	7642	7645	7648
7654	7656	7657	7660	7662	7677	7678	7683	7684	7688	7695	7701	7703	7705
7707	7711	7727	7729	7741	7744	7749	7750	7752	7762	7766	7769	7770	7780
7781	7785	7786	7788	7790	7794	7795	7798	7802	7803	7810	7818	7822	7825
7831	7835	7838	7840	7841	7843	7845	7858	7862	7868	7872	7884	7885	7886
7888	7900	7902	7903	7906	7918	7923	7924	7929	7932	7934	7936	7938	7942
7948	7954	7963	7968	7969	7973	7976	7984	7987	7989	7992	7994	7997	7998
8001	8004	8008	8009	8012	8019	8022	8026	8030	8036	8039	8045	8049	8050
8054	8059	8061	8062	8063	8064	8066	8073	8075	8084	8086	8087	8098	8099
8101	8102	8105	8110	8111	8112	8116	8121	—	—	—	—	—	—

Notes. Temporal and spectral data for these triggers are available in the BATSE 4B and Current Catalogs. The spectral peak energy (E_p) estimates of the above triggers and the rest of the 2130 BATSE Catalog GRBs are provided by Shahmoradi & Nemiroff (2010). The full conditional E_p probability density functions are available for download at https://sites.google.com/site/amshportal/research/aca/in-the-news/lgrb-world-model.

Download table as:  ASCIITypeset images: 1 2

2.2. Model Construction

The goal of the presented analysis is to derive a multivariate model that is capable of reproducing the observational data of the 1366 BATSE LGRBs. Examples of multivariate treatment of LGRB data are rare in GRB literature, with the most recent (and perhaps the only) example of such work presented by B10. Conversely, many authors have focused primarily on the univariate distribution of the spectral parameters, most importantly on the LF. A variety of univariate models have been proposed as the LGRB LF and fit to data by approximating the complex detector threshold as a step function (Schmidt 1999) or an efficiency grid (e.g., the four-interval efficiency modeling of Guetta et al. 2005) or by other approximation methods. A more accurate modeling of the LF, however, requires at least two LGRB observables incorporated in the model: the bolometric peak flux (P_bol) and the observed peak energy (E_p). The parameter E_p is required, since most gamma-ray detectors are photon counters, a quantity that depends on not only P_bol but also E_p of the burst. This leads to the requirement of using a bivariate distribution as the minimum acceptable model to begin with, for the purpose of constraining the LF. The choice of model can be almost anything (e.g., Kommers et al. 2000; Porciani & Madau 2001; Sethi & Bhargavi 2001; Schmidt 2009; Campisi et al. 2010; Wang & Dai 2011), since current theories of LGRB prompt emission do not set strong limits on the shape and range of the LF or any other LGRB spectral or temporal variables.

Here, the multivariate log-normal distribution is proposed as the simplest natural candidate model capable of describing data. The motivation behind this choice of model comes from the available observational data that closely resemble a joint multivariate log-normal distribution for the four most widely studied temporal and spectral parameters of LGRBs in the observer frame: P_bol, S_bol (bolometric fluence), E_p, T₉₀: since most LGRBs originate from moderate redshifts z ∼ 1–3, a fact known thanks to Swift satellite (e.g., B10; Racusin et al. 2011), the convolution of these observer-frame parameters with the redshift distribution results in negligible variation in the shape of the rest-frame joint distribution of the same LGRB parameters. Therefore, the redshift-convoluted four-dimensional (4D) rest-frame distribution can be well approximated as a linear translation from the observer-frame parameter space to the rest-frame parameter space, keeping the shape of the distribution almost intact. This implies that the joint distribution of the intrinsic LGRB variables: the isotropic peak luminosity (L_iso), the total isotropic emission (E_iso), the rest-frame time-integrated spectral peak energy (E_{p, z}), and the rest-frame duration (T_{90, z}) might be indeed well described by a multivariate log-normal distribution.

In general, models with higher nonzero moments than the log-normal model can also be considered for fitting, such as a multivariate skew-lognormal (e.g., Azzalini 1985) or variants of multivariate stable distributions (e.g., Press 1972). This, however, requires fitting for a higher number of free parameters, which is practically impossible given BATSE data with no available redshift information.

Following the discussion above, the process of LGRB observation can be therefore considered as a non-homogeneous Poisson process whose mean rate parameter—the cosmic LGRB differential rate, $\mathcal {R}_{{\rm cosmic}}$—is the product of the differential comoving LGRB rate density $\dot{\zeta }(z)$ with a p = 4D log-normal probability density function (pdf), $\mathcal {LN}$, of four LGRB variables: L_iso, E_iso, E_{p, z}, and T_{90, z}, with location vector μ and the scale (i.e., covariance) matrix Σ,

$$\begin{eqnarray} \mathcal{R}_{\rm cosmic} &=& \frac{\it dN}{{\it dL}_{\rm iso}\,{\it dE}_{\rm iso}\,{\it dE}_{p,z}\,{\it dT}_{90,z}\,{\it dz}} \\ &\propto & \mathcal {LN} (L_{{\rm iso}},E_{{\rm iso}},E_{p,z},T_{90,z}|\boldsymbol \mu,\Sigma ) \nonumber \\ &\times & \frac{\dot{\zeta }(z){{\it dV}/{\it dz}}}{(1+z)}, \nonumber \end{eqnarray} \tag{ 1 }$$

where the factor (1 + z) in the denominator accounts for cosmological time dilation and the comoving volume element per unit redshift, dV/dz, is

$$\begin{equation} \frac{\it dV}{\it dz} = \frac{C}{H_0}\frac{4\pi {D_L}^2(z)}{(1+z)^2 [\Omega _M(1+z)^3+\Omega _\Lambda ]^{1/2}}, \end{equation} \tag{ 2 }$$

with D_L being the luminosity distance,

$$\begin{equation} D_L(z)=\frac{C}{H_0}(1+z)\int ^{z}_{0}dz^{\prime } [(1+z^{\prime })^{3}\Omega _{M}+\Omega _{\Lambda } ]^{-1/2}, \end{equation} \tag{ 3 }$$

assuming a flat ΛCDM cosmology, with parameters set to h = 0.70, Ω_M = 0.27, and $\Omega _\Lambda =0.73$ (Jarosik et al. 2011). Here, C and H₀ = 100 h (Km s⁻¹ MPc⁻¹) stand for the speed of light and the Hubble constant, respectively.

The 4D log-normal distribution of Equation (1), $\mathcal {LN}$, has an intimate connection to multivariate Gaussian distribution in the logarithmic space of LGRB observable parameters (cf. Appendix D).

One could generalize the LGRB rate of Equation (1) to incorporate a redshift evolution of the LGRB variables in the form of μ_i(z) = μ_{0, i} + α_ilog (1 + z), (i = 1, ..., 4), where α has to be constrained by observational data. This is, however, impractical for BATSE data due to unknown redshifts, as the fitting results in degenerate values for α. Nevertheless, the multivariate analysis of Swift LGRBs presented by B10 strongly rejects the possibility of redshift evolution of the LF, a fact that further legitimizes the absence of redshift–luminosity evolution in Equation (1).

As for the comoving rate density $\dot{\zeta }(z)$, it is assumed that LGRBs trace the star formation rate (SFR) in the form of a piecewise power-law function of Hopkins & Beacom (2006, hereafter HB06):

$$\begin{equation} \dot{\zeta }(z) = \frac{\it dN}{\it dz} \propto \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(1+z)^{\gamma _0} & z<z_0 \\ (1+z)^{\gamma _1} & z_0<z<z_1 \\ (1+z)^{\gamma _2} & z>z_1, \\ \end{array}\right. \end{equation} \tag{ 4 }$$

with parameters (z₀, z₁, γ₀, γ₁, γ₂) set to best-fit values (0.97, 4.5, 3.4, −0.3, −7.8) of HB06, and also to the best values (0.993, 3.8, 3.3, 0.055, −4.46) of an updated SFR fit by Li (2008). Alternatively, the bias-corrected redshift distribution of LGRBs derived from Swift data (B10) with best-fit parameter values (0.97, 4.00, 3.14, 1.36, −2.92) can be employed as $\dot{\zeta }(z)$. This parameter set is consistent with an LGRB rate scenario tracing metallicity-corrected SFR with a cutoff Z/Z_☉ ∼ 0.2–0.5 (Figure 10 and Equation (8) in B10; Li 2008). The hypothesis of LGRB rate evolving with cosmic metallicity is both predicted by the Collapsar model of LGRBs (e.g., Woosley & Heger 2006) and supported by observations of LGRB host galaxies (e.g., Stanek et al. 2006; Levesque et al. 2010a), although the metallicity–rate connection and the presence of a sharp metallicity cutoff have been challenged by few recent host galaxy observations (e.g., Levesque et al. 2010c, 2010b) and possible unknown observational biases (e.g., Levesque 2012).

The cosmic LGRB rate, $\mathcal {R}_{{\rm cosmic}}$, in Equation (1), although quantified correctly, does not represent the observed rate ($\mathcal {R}_{{\rm obs}}$) of LGRBs detected by BATSE LADs, unless convolved with an accurate model of BATSE trigger efficiency, η, as a function of the burst redshift and rest-frame parameters, discussed in Appendix B,

$$\begin{equation} \mathcal {R}_{{\rm obs}}=\eta (L_{{\rm iso}},E_{p,z},T_{90,z},z)\times \mathcal {R}_{{\rm cosmic}}. \end{equation} \tag{ 5 }$$

2.3. Model Fitting

Now, with a statistical model for the observed rate of LGRBs at hand (i.e., Equation (5)), the best-fit parameters can be obtained by the method of maximum likelihood. This is done by maximizing the likelihood function of the model, given the observational data, using a variant of the Metropolis–Hastings Markov chain Monte Carlo (MCMC) algorithm discussed in detail in Appendix C. As mentioned before, the fitting is performed for three redshift-distribution scenarios of HB06, B10, and Li (2008).

It is also known that the T₉₀ of LGRBs are potentially subject to estimation biases. To ensure that the reported T₉₀ of BATSE LGRBs do not bias the fitting results for the rest of the parameters, model fitting was also performed by considering only three spectral variables of BATSE LGRBs: the bolometric 1 s peak flux (P_bol), bolometric fluence (S_bol), and observed peak energy (E_p), excluding duration (T_{90, z}) variable from the model, thus reducing the dimension of the model by one. Only after the fitting was performed did it become clear that the inclusion of the T₉₀ of BATSE LGRBs in the fitting does not significantly affect the resulting best-fit parameters of the model. Therefore, only results from the full model fitting are presented here, as in Table 2.

Table 2. Mean Best-fit Parameters of LGRB World Model, for the Three Redshift Distribution Scenarios

Parameter	HB06	Li (2008)	B10
Redshift parameters (Equation (4))
z0	0.97	0.993	0.97
z1	4.5	3.8	4.00
γ0	3.4	3.3	3.14
γ1	−0.3	0.0549	1.36
γ2	−7.8	−4.46	−2.92
Location parameters
log (Liso)	51.35 ± 0.20	51.50 ± 0.19	51.73 ± 0.19
log (Eiso)	51.82 ± 0.20	51.94 ± 0.20	52.03 ± 0.21
log (Ep, z)	2.43 ± 0.05	2.47 ± 0.05	2.54 ± 0.05
log (T90, z)	0.99 ± 0.03	0.96 ± 0.03	0.80 ± 0.03
Scale parameters
$\log (\sigma _{L_{{\rm iso}}})$	−0.22 ± 0.06	−0.23 ± 0.06	−0.20 ± 0.05
$\log (\sigma _{E_{{\rm iso}}})$	−0.07 ± 0.03	−0.07 ± 0.03	−0.04 ± 0.03
$\log (\sigma _{E_{p,z}})$	−0.44 ± 0.02	−0.44 ± 0.02	−0.44 ± 0.02
$\log (\sigma _{T_{90,z}})$	−0.38 ± 0.01	−0.39 ± 0.01	−0.40 ± 0.01
Correlation coefficients
$\rho _{L_{{\rm iso}}-E_{{\rm iso}}}$	0.93 ± 0.01	0.94 ± 0.01	0.96 ± 0.01
$\rho _{L_{{\rm iso}}-E_{p,z}}$	0.47 ± 0.07	0.45 ± 0.07	0.44 ± 0.08
$\rho _{L_{{\rm iso}}-T_{90,z}}$	0.52 ± 0.08	0.59 ± 0.09	0.75 ± 0.07
$\rho _{E_{{\rm iso}}-E_{p,z}}$	0.58 ± 0.04	0.58 ± 0.04	0.59 ± 0.04
$\rho _{E_{{\rm iso}}-T_{90,z}}$	0.63 ± 0.05	0.66 ± 0.05	0.74 ± 0.04
$\rho _{E_{p,z}-T_{90,z}}$	0.34 ± 0.04	0.37 ± 0.04	0.50 ± 0.04
BATSE LGRB detection efficiency (Equation (A5))
μthresh	−0.44 ± 0.02	−0.45 ± 0.02	−0.44 ± 0.02
log (σthresh)	−0.88 ± 0.05	−0.90 ± 0.05	−0.88 ± 0.05

Notes. The full Markov chain sampling of the above parameters from the 16-dimensional parameter space of the likelihood function (Appendix C) is available for download at https://sites.google.com/site/amshportal/research/aca/in-the-news/lgrb-world-model for each of the three redshift distributions.

Download table as:  ASCIITypeset image

Due to lack of redshift information for BATSE GRBs, the resulting parameters of the model exhibit strong covariations with each other. This is illustrated in the example correlation matrix of the LGRB world model in Table 3. All location parameters appear to correlate strongly positively with each other, and so do the scale parameters. The location parameters however negatively correlate with the scale parameters, meaning that an increase in the average values of the rest-frame parameters reduces the half-width of the corresponding distributions of the variables. In general, it is also observed that the correlations among the four variables weaken with increasing the location parameters. An exception to this is the correlation of T_{90, z} with L_iso and E_iso which tends to increase with location parameters. Since an excess in the cosmic rates of LGRBs at high redshifts generally results in an increase in the values of location parameters, it can be said that "given BATSE LGRBs data, a higher rate for LGRBs at distant universe generally implies weaker E_iso–E_{p, z} and L_iso–E_{p, z} correlations and stronger covariation of T_{90, z} with the three other parameters."

Table 3. Correlation Matrix of the Parameters of the LGRB World Model, for the Median Case of an LGRB Cosmic Rate Tracing Star Formation Rate in Li (2008)

Parameter	log (Liso)	log (Eiso)	log (Ep, z)	log (T90, z)	$\log (\sigma _{L_{{\rm iso}}})$	$\log (\sigma _{E_{{\rm iso}}})$	$\log (\sigma _{E_{p,z}})$	$\log (\sigma _{T_{90,z}})$	$\rho _{L_{{\rm iso}}-E_{{\rm iso}}}$	$\rho _{L_{{\rm iso}}-E_{p,z}}$	$\rho _{L_{{\rm iso}}-T_{90,z}}$	$\rho _{E_{{\rm iso}}-E_{p,z}}$	$\rho _{E_{{\rm iso}}-T_{90,z}}$	$\rho _{E_{p,z}-T_{90,z}}$	μthresh	log (σthresh)
log (Liso)	1.00	0.99	0.90	0.34	−0.91	−0.86	−0.59	−0.14	−0.10	−0.45	0.51	−0.52	0.45	0.05	−0.68	−0.44
log (Eiso)		1.00	0.92	0.42	−0.91	−0.90	−0.62	−0.15	−0.20	−0.54	0.45	−0.56	0.40	0.00	−0.67	−0.43
log (Ep, z)			1.00	0.38	−0.82	−0.84	−0.79	−0.15	−0.26	−0.77	0.40	−0.77	0.38	0.07	−0.61	−0.40
log (T90, z)				1.00	−0.37	−0.52	−0.32	−0.16	−0.52	−0.32	−0.48	−0.33	−0.50	−0.56	−0.17	−0.09
$\log (\sigma _{L_{{\rm iso}}})$					1.00	0.94	0.59	0.12	0.14	0.53	−0.53	0.57	−0.46	0.03	0.50	0.30
$\log (\sigma _{E_{{\rm iso}}})$						1.00	0.67	0.17	0.41	0.63	−0.33	0.66	−0.32	0.14	0.46	0.28
$\log (\sigma _{E_{p,z}})$							1.00	0.12	0.36	0.84	−0.22	0.84	−0.25	−0.06	0.37	0.24
$\log (\sigma _{T_{90,z}})$								1.00	−0.03	0.10	−0.09	0.11	0.00	0.03	0.11	0.07
$\rho _{L_{{\rm iso}}-E_{{\rm iso}}}$									1.00	0.47	0.43	0.41	0.22	0.26	−0.01	−0.01
$\rho _{L_{{\rm iso}}-E_{p,z}}$										1.00	−0.15	0.97	−0.21	−0.05	0.29	0.19
$\rho _{L_{{\rm iso}}-T_{90,z}}$											1.00	−0.18	0.95	0.57	−0.32	−0.20
$\rho _{E_{{\rm iso}}-E_{p,z}}$												1.00	−0.22	0.06	0.29	0.18
$\rho _{E_{{\rm iso}}-T_{90,z}}$													1.00	0.64	−0.26	−0.17
$\rho _{E_{p,z}-T_{90,z}}$														1.00	−0.09	−0.07
μthresh															1.00	0.79
log (σthresh)																1.00

Download table as:  ASCIITypeset image

2.4. Goodness-of-fit Tests

In any statistical fitting problem, perhaps more important than the model construction is to provide tests showing how good the model fit is to input data. For many univariate studies of the GRB LF, this is done by employing well-established statistical tests such as the Kolmogorov–Smirnov (K-S; e.g., Kolmogoroff 1941; Smirnov 1948) or Pearson's χ² (e.g., Fisher 1924) tests. In the case of multivariate studies (e.g., B10), a combination of visual inspection of the fitting results, K-S test on the marginal, and bivariate distributions and variants of χ² (e.g., likelihood ratio) tests have been used.

In general, univariate tests on the marginal distributions of multivariate fits provide only necessary—but not sufficient—evidence for a good multivariate fit. Alternatively one could assess the similarity by using nonparametric multivariate GoF tests. Such tests, although existing, have been rarely discussed and treated in statistics due to difficulties in the interpretation of the test statistic (e.g., Peacock 1983; Press et al. 1992; Justel et al. 1997). Ideally, one can always use the Pearson's χ² GoF test for any multivariate distribution. However, for the special case of BATSE 1366 LGRBs, one would need an observed sample consisting of N ≫ 1366 observations to avoid serious instabilities that occur in χ² tests due to small sample sizes (e.g., Cochran 1954).

To ensure a good fit to the observational data in all—and not only univariate—levels of the multivariate structure of data, an assessment of similarity can be obtained by scanning and comparing the model and data along their principal axes, in addition to univariate tests on the marginal distributions. Although statistically not a sufficient condition for the multivariate similarity of the model prediction to data, this can provide strong evidence in favor of a good fit, at a much higher confidence than tests performed only on the marginal distributions.

Following from above, Figure 2 presents the model predictions for marginal distributions of the four LGRB variables in the observer frame. The K-S test probabilities for the similarity of the model predictions to the marginal distributions of BATSE LGRB variables are also reported on the top right of each plot. All three redshift-distribution parameters of HB06, Li (2008), and B10 (Equation (4)) result in relatively similar fits to BATSE data in the observer frame. Thus, for brevity, only plots for one representative (median) redshift distribution (i.e., Li 2008) are presented.

Zoom In Zoom Out Reset image size

As the second level of GoF tests, the joint bivariate model predictions are compared to BATSE LGRB data, presented in Figures 3–5. This method of scanning model and data along the principal axes of the joint bivariate distributions can be generalized to trivariate and quadruvariate joint distributions. For brevity, however, only the bivariate tests are presented.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is observed in the plots of Figures 2–5 that the model provides excellent fit to data, within the uncertainties caused by random Poisson fluctuations in the BATSE LGRB observed rate. These random fluctuations in BATSE detections are encompassed in each graph by the green dashed lines that represent the 90% confidence intervals (CI) on BATSE LGRB detections (blue solid lines), derived by repeated sampling from the model.

Unfortunately, the same methods for a comparison of data and model cannot be applied in the LGRB rest frame, due to lack of redshift information for the BATSE sample of LGRBs. Nevertheless, a comparison of the model with observational data of other instruments—with measured redshifts—can provide clues on the underlying joint distribution of LGRB temporal and spectral variables in the rest frame compared to LGRB detections of different gamma-ray instruments, as will be done in the following sections.

3.1. LGRB Luminosity Function and log(N)–log(P) Diagram

The log (N)–log (P) diagram of GRBs has been subject of numerous studies in the BATSE era, primarily for the purpose of finding signatures of cosmological (versus galactic) origins in the LGRB rate. The cosmological origin of LGRBs is now well established. Nevertheless, the log (N)–log (P) diagram can still provide useful information for future gamma-ray experiments.

Figure 2 (bottom) depicts the prediction of the LGRB world model for the traditional log (N)–log (P) diagram for 1 s peak photon flux in the BATSE nominal detection energy range 50–300 keV, both for the differential (left panel) and the cumulative (right panel) LGRB rate. For all three LGRB cosmic rates considered in this work—as in Table 2—the differential log (N)–log (P) diagram shows a peak in the rate at P_50–300 ∼ 0.1 photons s⁻¹ cm⁻². Such a peak in the LGRB rate results in a relative flattening at the dim end of the cumulative log (N)–log (P) diagram, as compared to its bright end. This observation has already been reported by B10 for the Swift sample of LGRBs, although an entirely different LF—a broken power-law LF—was used by B10 in their multivariate LGRB world model (cf. Equation (2) in B10).

On the other hand, the peak in the observer-frame LGRB rate translates to three relatively different (at ∼1σ level) peaks in the LF of LGRBs. In general, it is observed that the peak of the LF (i.e., the average 1 s peak luminosity of LGRBs) increases with increasing cosmic rate of LGRBs at high redshift. This effect is well depicted in the top left panel of Figure 6 for the three LGRB cosmic rates considered: HB06, Li (2008), and B10.

Zoom In Zoom Out Reset image size

Compared to the predictions of B10's LGRB world model (cf. bottom plot of Figure 6 in B10), the log-normal model suggests a lower peak for BATSE LGRB LF (∼52.3 here versus ∼52.7 in B10) for the same redshift distribution of LGRBs. Averaging over the three redshift distributions considered, the model predicts a dynamic 3σ range of observer-frame brightness log (P_bol(erg s⁻¹ cm⁻²)) ∈ [ − 7.11 ± 2.66] corresponding to P_bol(erg s⁻¹ cm⁻²) ∈ [1.70 × 10⁻¹⁰, 3.58 × 10⁻⁵] for LGRBs. This translates to an average dynamic 3σ range—in the rest frame—of log (L_iso(erg s⁻¹)) ∈ [51.53 ± 1.99] corresponding to L_iso(erg s⁻¹) ∈ [3.46 × 10⁴⁹, 3.38 × 10⁵³].

3.2. Isotropic Emission and Peak Energy Distributions

3.3. Duration Distribution

3.4. Temporal and Spectral Correlations

Ever since the launch of the Swift gamma-ray detector satellite, there has been a flurry of reports on the discovery of strong and significant correlations among the spectral parameters of LGRBs, most prominently among the rest-frame spectral peak energy and the total isotropic emission or peak luminosity of LGRBs (e.g., Amati et al. 2008; Ghirlanda et al. 2008). Despite the lack of measured redshifts for BATSE GRBs, signatures of such correlations had been found by earlier works in the BATSE era through careful analysis of observer-frame spectral properties of LGRBs (e.g., Lloyd et al. 2000). Nevertheless, the strength and significance of these correlations were undermined by analyzing larger samples of the BATSE catalog of GRBs (e.g., Nakar & Piran 2005; Band & Preece 2005; Shahmoradi & Nemiroff 2009, 2011) or the Swift sample of GRBs (e.g., Butler et al. 2007, 2009; B10), all arguing that the sample of bursts used to construct the claimed spectral relations is representative of only bright-soft LGRBs. These arguments have been responded to by others (e.g., Ghirlanda et al. 2005, 2008; Nava et al. 2008; see Shahmoradi & Nemiroff 2011 for a complete history of the debate).

A multivariate analysis of the BATSE LGRB data that carefully eliminates potential biases at the detection process can therefore greatly help us understand the strength and significance of the reported correlations. Figure 7 shows the predictions of the LGRB world model for the two widely discussed spectral correlations: the L_iso–E_{p, z} (the Yonetoku) and E_iso–E_{p, z} (the Amati) relations. As indicated by the model, a large fraction of BATSE LGRBs (and a much larger fraction of the entire LGRB population) on the dim-hard regions of the plots appear to be underrepresented by the sample of LGRBs used for the construction of these relations. Nevertheless, given the three redshift distributions considered for the model, a relatively strong (Pearson's correlation coefficient $\rho _{E_{{\rm iso}}-E_{p,z}}=0.58\pm 0.04$) and highly significant (>14σ) correlation is predicted between the E_iso and E_{p, z} of LGRBs. The slope of the two relations suggested by the model also differs significantly from the original reports of the relations (Schaefer 2007; Ghirlanda et al. 2008). Generally, in regression modeling, it is assumed that there are two independent and dependent variables. In the case of the proposed relations, none of the variables is known to depend theoretically on the other. Therefore, due to the large statistically unexplained variances of the two LGRB variables, different regression methods, such as ordinary least squares, OLS(Y|X) and OLS(X|Y), result in entirely different slopes for the relations. The best-fit power-law relations and the conditional variances of the regress and given the regressor can be easily obtained from the parameters of the multivariate log-normal model in Table 2 (e.g., Section 2.11 in Kutner et al. 2004).

Zoom In Zoom Out Reset image size

A partial correlation analysis of the two L_iso–E_{p, z} and E_iso–E_{p, z} relations (Figure 8) reveals that the moderate correlation of the isotropic peak luminosity with the time-integrated peak energy is entirely due to the strong association of L_iso with the total isotropic emission from the burst. As seen in the top left plot of Figure 8 there is indeed a negative correlation between the L_iso and E_{p, z} of GRBs for a fixed isotropic emission E_iso and burst duration T_{90, z}. Conversely, the model indicates a highly significant correlation of E_iso with the time-integrated E_{p, z}, even after elimination of the effects of L_iso and T_{90, z} from the E_iso–E_{p, z} relation.

Zoom In Zoom Out Reset image size

As observed in Table 2, the model predicts positive correlation among all four LGRB variables. The isotropic peak luminosity (L_iso) and the total isotropic emission (E_iso) appear to be strong indicators of each other reciprocally. Surprisingly, it is also observed that the rest-frame duration T_{90, z} of LGRBs strongly correlates with both E_iso and L_iso. The existence of a possible positive correlation between the isotropic emission and the duration of LGRBs has been implied by the analysis of Swift LGRBs (B10), though only weakly present there. Such correlations can be enlightening for the early studies of time-dilation signatures in BATSE GRBs. A positive duration–brightness correlation is also opposite to—but not necessarily in contradiction with—the negative duration–brightness correlation in pulses of individual GRBs (e.g., Fenimore et al. 1995; Nemiroff 2000; Ramirez-Ruiz & Fenimore 2000). Combination of the two correlations implies that the number of pulses in individual LGRBs should be positively correlated with the peak luminosity (or equivalently, the total isotropic emission) of the bursts. This is indeed in qualitative agreement with the observed inequality relation between the isotropic peak luminosity and the number of pulses in Swift LGRBs (Figure 6 in Schaefer 2007). The strength of the correlations found encourages the search for the underlying physical mechanism that could give rise to these relations. This is, however, beyond the scope of this paper (cf. Rees & Mészáros 2005; Ryde et al. 2006; Thompson et al. 2007; Giannios 2012; Dado & Dar 2012 for example discussions).

It is also worth mentioning that the T₉₀ duration of BATSE LGRBs strongly correlates with the bolometric fluence to bolometric 1 s peak flux ratio (FPR), with Pearson's correlation coefficient ρ_FPR–T90 ≈ 0.67. A comparison of BATSE data with the predictions of the model for the bivariate distribution of FPR and T₉₀ is given in Figure 9 (left panel). Interestingly, the model predicts the same correlation strength of ρ_FPR–T90 ≈ 0.67 for the entire LGRB population, implying that the detection process does not bias the FPR–T₉₀ relation. Such strong correlation indicates an underlying intrinsic interrelation between the three variables: P_bol, S_bol, and T₉₀, also among their corresponding rest-frame counterparts. In fact, B10 use a variant of this trivariate correlation to define an effective peak flux in terms of fluence and duration, discarding the traditional definition of peak flux as the peak photon counts in 1 s time interval.

Zoom In Zoom Out Reset image size

The primary goal of the presented analysis was to model and constrain the LF, temporal and spectral correlations, and energetics of long-duration class of GRBs by exploiting the wealth of information that has been buried and untouched in the BATSE GRB catalog to this date. Below is a summary of steps taken to construct the LGRB world model, detailed in Section 2.

• 1.  
  A sample of 1366 LGRBs is carefully selected by using the fuzzy C-means clustering algorithm based on temporal and spectral parameters and inspection of the individual light curves of BATSE-catalog GRBs (Figure 1; Section 2.1).
• 2.  
  It is proposed that the BATSE LGRB data might be very well consistent with being drawn from a multivariate log-normal population of LGRBs in four rest-frame LGRB variables: the bolometric isotropic 1 s peak luminosity (L_iso), the bolometric isotropic emission (E_iso), the spectral peak energy (E_{p, z}), and the duration (T_{90, z}). Therefore, the observed joint distribution of the four LGRB variables: the bolometric 1 s peak flux (P_bol), the bolometric fluence (S_bol), the observed spectral peak energy (E_p) and the observed duration (T₉₀) results from the convolution of the rest-frame multivariate log-normal population with the cosmic rate (i.e., the redshift distribution) of LGRBs, truncated by the complex LGRB trigger threshold of BATSE LADs, as illustrated in Section 2.2, Equations (1)–(5). A prescription for modeling BATSE LAD detection efficiency is given in Appendix B.
• 3.  
  The LGRB model (Equation (1)) is fit to BATSE data by maximizing the likelihood function of the model (Section 2.3 and Appendix C; Equation (C2)). In order to derive the best-fit parameters of the model and their corresponding uncertainties, an adaptive Metropolis–Hastings MCMC (AMH-MCMC) algorithm is set up to efficiently sample from the 16D likelihood function. The best-fit parameters are obtained for three LGRB cosmic rates: SFR of HB06, SFR of Li (2008), and the predicted LGRB redshift distribution of B10 which is consistent with the LGRB rate tracing cosmic metallicity with a cutoff Z/Z_☉ ∼ 0.2–0.5.
• 4.  
  To ensure the model provides adequate fit to observational data, multivariate GoF tests are presented (Section 2.4 and Figures 2–5).

Summarized below are the principal conclusions drawn from the analysis based on the proposed LGRB world model.

• 1.  
  Energetics. It is expected that the peak brightness distribution of LGRBs has the effective range of log (P_bol(erg s⁻¹ cm⁻²)) ∈ [ − 7.11 ± 2.66] corresponding to P_bol(erg s⁻¹ cm⁻²) ∈ [1.70 × 10⁻¹⁰, 3.58 × 10⁻⁵]. This translates to a dynamic 3σ range—in the rest frame—of log (L_iso(erg s⁻¹)) ∈ [51.53 ± 1.99], corresponding to L_iso(erg s⁻¹) ∈ [3.46 × 10⁴⁹, 3.38 × 10⁵³]. In addition, a turnover is predicted in the differential log (N)–log (P) diagram of LGRBs at P_50–300 ∼ 0.1 photons s⁻¹ cm⁻² in the BATSE nominal detection energy range (50–300 keV). This is consistent with and further extends the apparent flattening in the cumulative log (N)–log (P) diagram of Swift LGRBs reported recently by B10.As for the bolometric fluence and the total isotropic emission distributions, a range of log (S_bol(erg cm⁻²)) ∈ [ − 6.16 ± 3.01] corresponding to S_bol(erg cm⁻²) ∈ [6.82 × 10⁻¹⁰, 7.01 × 10⁻⁴] is indicated. This translates to an average dynamic 3σ range—in the rest frame—of log (E_iso (erg)) ∈ [51.93 ± 2.71] corresponding to E_iso (erg) ∈ [1.66 × 10⁴⁹, 4.46 × 10⁵⁴] (Sections 3.2 and 3.1; Table 2).
• 2.  
  Durations and spectral peak energies. The rest-frame spectral peak energies (E_{p, z}) of LGRBs are likely well described by a log-normal distribution with an average 3σ range of log (E_{p, z} (keV)) ∈ [2.48 ± 1.12] corresponding to E_{p, z} (keV) ∈ [23, 4006] with peak LGRB rate at E_{p, z} ∼ 300 (keV). This translates to an effective observer-frame peak energy range of log (E_p (keV)) ∈ [1.93 ± 1.22] corresponding to E_p (keV) ∈ [5, 1427] with peak LGRB rate at E_p ∼ 85 keV. It is also observed that the observer-frame T₉₀ duration of LGRBs peaks at T₉₀ ∼ 30 s with a 3σ range of T₉₀ (s) ∈ [1.4, 620]. This translates to an average 3σ range of rest-frame log (T_{90, z} (s)) ∈ [0.92 ± 1.24] corresponding to T_{90, z} (s) ∈ [0.47, 145] with a peak rate at T_{90, z} ∼ 10 s (Section 3.3; Table 2).Recently, Bromberg et al. (2012) proposed the apparent flatness in the duration distribution of BATSE LGRBs—when plotted in the form of dN/dT₉₀ instead of dN/dlog (T₉₀)—as the first direct evidence of the Collapsar model of LGRBs. The results of presented analysis are inconsistent with a flat T₉₀ distribution of LGRBs at short durations (Figure 2, center right panel and Figure 9, right panel). The observed flat T₉₀ distribution of LGRBs at short durations can be explained away in terms of the skewed nature of the log-normal distribution subject to sample incompleteness. It is therefore expected that a significantly larger sample of LGRBs that will be detected by future gamma-ray satellites will smear out the apparent flatness at the short tail of the duration distribution of LGRBs. A similar flat distribution is also observed for SGRBs at very short durations (Figure 9, right panel) which might be hard to reconcile with the Collapsar interpretation of the observed flatness in the LGRB T₉₀ distribution, proposed by Bromberg et al. (2012).
• 3.  
  Temporal and spectral correlations. All four LGRB variables: L_iso, E_iso, E_{p, z}, and T_{90, z} appear to be either moderately or strongly positively correlated with each other. In particular, a relatively strong and "broad" but highly significant correlation strength (Pearson's correlation coefficient $\rho _{E_{{\rm iso}}-E_{p,z}}=0.58\pm 0.04$) is predicted between E_iso and E_{p, z} of long-duration class of GRBs. Surprisingly, T_{90, z} appears to evolve with L_iso and E_iso such that brighter bursts generally tend to have longer durations (Section 3.4; Table 2). This prediction of the model together with the previously reported negative correlation of the brightness and the duration of individual pulses in LGRBs (e.g., Fenimore et al. 1995; Nemiroff 2000; Ramirez-Ruiz & Fenimore 2000) might possibly indicate that intrinsically brighter LGRBs contain, on average, a higher number of pulses.There is a slight chance that a small fraction (<50) of BATSE LGRBs were misclassified as SGRBs by the automated pattern recognition methods exploited in this analysis (see Figure 3, center right panel). If true, it will most likely affect (if significant at all) the constraints derived on the LF of LGRBs and the correlation of E_{p, z} with L_iso.
• 4.  
  Redshift distribution. The lack redshift information for the BATSE GRBs strongly limits the prediction power of the presented analysis for the cosmic rate of LGRBs. Nevertheless, based on the Markov chain sampling of the likelihood function for the three LGRB redshift distributions considered here (Section 2.2 and Figure 10), it is observed that BATSE data potentially, but not necessarily, favor an LGRB rate consistent with cosmic metallicity evolution with a cutoff Z/Z_☉ ∼ 0.2–0.5 (cf. B10), with no luminosity–redshift evolution.Assuming that LGRBs track SFR, only a tiny fraction (i.e., ∼2–3) of 1366 BATSE LGRBs are expected to have originated from high redshifts (z ≳ 5). In the case of an LGRB rate tracing cosmic metallicity evolution (e.g., B10), the fraction increases by one order of magnitude to ∼2%, corresponding to ∼27 bursts out of 1366 BATSE LGRBs. For comparison, the expected fraction of Swift and EXIST LGRBs with z ≳ 5 are ∼6% and ∼7% (B10). The discrepancy is well explained by the fact that both Swift and EXIST are more sensitive to long-soft bursts—characteristic of high-redshift LGRBs—due to their lower gamma-ray trigger energy window, compared to BATSE DADs (cf. Gehrels et al. 2004; Band et al. 2008; Grindlay & the EXIST Team 2009).

Zoom In Zoom Out Reset image size

Although fitting is performed for the rest-frame variables, it is notable that the overall shape of the resulting observer-frame distribution of the variables also resembles a multivariate log-normal (see Figures 2–7). In other words, the redshift convolution of the rest-frame population distribution approximately acts as a linear transformation from the rest frame of LGRBs to the observer frame. This is primarily due to the narrow redshift distribution of LGRBs—as compared to the width of the LGRB rest-frame temporal and spectral distributions—with almost 90% of the population originating from intermediate redshifts, z ∈ [1, 4.3]. Balazs et al. (2003) provide an elegant discussion on the potential effects of redshift convolution on the observed distribution of LGRB durations and spectral parameters.

As implied by the model, there is no evidence for a significant population of bright-hard LGRBs that could have been missed in the BATSE catalog of GRBs. Conversely, a large population of low luminosity with moderate-to-low spectral peak energies seem to have gone undetected by BATSE LAD detectors. It should be emphasized that the apparent lack of very bright-soft LGRBs has a true physical origin according to the analysis presented here and is not an artifact of the detection process or spectral fitting models (e.g., the Band model, CPL or SBPL models) used by GRB researchers. Whether the X-Ray Flashes (XRF), X-ray-rich, and the sub-luminous GRBs (e.g., Strohmayer et al. 1998; Kippen et al. 2003; Tikhomirova et al. 2006) can be incorporated into a unified class of events described by a single model remains an open question in this work. At present, these events can be either considered as a separate class of cosmological events or as the soft-dim tail of the LGRB world model that have been mostly out of BATSE detection range and missed (see Figure 3 of Kippen et al. 2003 for a comparison with the predictions of the LGRB world model here in Figures 2–5). A definite answer to this question requires knowledge of the true rate of sub-luminous bursts and XRFs based on the observed rates of these events convolved with complex detection thresholds of different instruments used for observations.

B10 have presented an elaborate multivariate analysis of Swift LGRBs. While providing reasonable fit to Swift data, the model of B10 is primarily aimed at the discovery of the potential sub-luminous events that mostly go undetected by gamma-ray detectors. Such a model, capable of accounting for a large population of undetected sub-luminous bursts, is proposed at the cost of throwing away parts of information stored in the spectral parameters of LGRBs in the analysis of B10 (cf. Ghirlanda et al. 2012). Nevertheless, the presented analysis indicates that the apparent correlation of the isotropic peak luminosity (L_iso) with the time-integrated spectral peak energy (E_{p, z}) of LGRBs is peripheral to the more fundamental relation between the total isotropic emission (E_iso) and E_{p, z}, and that the relation can be created by defining an effective luminosity based on the two GRB variables E_iso and duration (e.g., T_{90, z}) (Figures 8 and 9, left panel). In light of the analysis presented by Ghirlanda et al. (2012), it can be therefore suggested that a new definition of luminosity based on the E_iso and T_{90, z} of GRBs drawn from the LGRB world model of B10 will alleviate the apparent discrepancy between the observed L_iso–E_{p, z} relation of LGRBs and the predicted relation from the LGRB model of B10. It is also expected that a better definition of peak luminosity that is not limited to a specific timescale in the observer frame of the bursts would result in an improvement in the correlations of the time-integrated spectral peak energy (E_{p, z}) and the isotropic emission (E_iso) with the luminosity variable (L_iso). Given the above arguments, the presence of L_iso as an independent variable in the LGRB world model—which was unfortunate due to the dependence of the BATSE trigger algorithm on the GRB's peak luminosity—might be viewed as overfitting and unnecessary.

The proposed multivariate log-normal model while requiring minimal free parameters compared to any other statistical model considered in GRB literature to date provides an accurate comprehensive description of the largest catalog of long-duration GRBs, serving as a powerful probe to explore the population properties of a large fraction of LGRBs that are missed in spectral analyses due to low-quality data or the lack of measured redshift, or simply go undetected due to the instrument's gamma-ray detector threshold. Data from future gamma-ray experiments will enable us to further confirm, improve, or invalidate predictions of the presented model.

The author owes special thanks to Robert J. Nemiroff (NASA GSFC Astrophysicist and professor of Physics at Michigan Tech University) for his generous support and help through numerous discussions over more than a year, from the early stages of developing the basic ideas behind this work to the core of the analysis and final steps of manuscript preparation.

I thank William H. Press (professor of Computer Science and Integrative Biology at The University of Texas at Austin), Lars Koesterke (at Texas Advanced Computing Center), and Mehdi Mortazavi (at MTU) for helpful comments on some statistical and computational aspects of this work. I am very grateful to Swadesh Mahajan (professor at Institute for Fusion Studies) and Richard Hazeltine (professor and chair at the Department of Physics, The University of Texas at Austin) for their generous support in the final stages of preparing this manuscript. I also thank Giancarlo Ghirlanda (at INAF-Osservatorio Astronomico di Brera) for his timely feedback and comment on the GRB terminology used in this manuscript. Greatly appreciated were helpful comments and detailed criticisms from the anonymous referee of this manuscript.

This work would have not been accomplished without the vast time and effort spent by many scientists and engineers who designed, built, and launched the Compton Gamma-Ray Observatory and were involved in the analysis of GRB data from BATSE Large Area Detectors.

It is well known that the traditional classification of GRBs based on a sharp cutoff in the observed duration (T₉₀) distribution of GRBs—usually set at T₉₀ = 2 s[50–300 keV]—is insufficient and can be misleading close to the sharply defined border. The apparent long-soft-bright to short-hard-dim trend observed in the prompt-emission properties of BATSE GRBs (e.g., Figure 1, top panel; see also Figure 8 of Shahmoradi & Nemiroff 2010) necessitates the use of a rigorous classification scheme based on all available spectral properties of GRBs, in addition to duration.

Despite a rich literature on the classification methodologies for GRBs (e.g., Hakkila et al. 2000a, 2003, and references therein), the choice of a classification method to separate the BATSE catalog of GRBs into two subgroups of long and short durations with minimal misclassifications remains a difficult task. This is primarily due to the significant overlap (or similarity) in all (or some) spectral and temporal properties of the two classes of GRBs, in addition to the heterogeneity of objective functions that might differ considerably from one classification algorithm to another.

Here, fuzzy (soft) clustering algorithms are preferred over hard clustering methods, since they provide a probability of the event belonging to each specific subgroup, in contrast to hard classifications that return only binary probabilities of either 0 or 1. The choice of fuzzy algorithm greatly facilitates identification of bursts that might have been potentially misclassified (Section 2.1).

Investigation of different fuzzy algorithms available in the literature leads us to two prominent candidates: the SAND method of Rousseeuw et al. (1996) and the fuzzy C-means discussed by Dunn (1973) and Bezdek (1981). While fuzzy C-means is especially useful for cases where subgroups are known to be approximately symmetric, the SAND algorithm is superior to C-means for its lack of sensitivity to different subgroup sizes, orientations, and asymmetries. Nevertheless, the presence of a handful of SGRs in the BATSE catalog—with spectral properties comparable to that of LGRBs—results in relatively poor classification by the SAND method as compared to C-means. Besides the choice of algorithm, the GRB variables—by which the classification is done—are selected such that the resulting relative sizes of the two SGRB and LGRB populations correspond to those found by Shahmoradi & Nemiroff (2010) through a different approach that the author believes to be less prone to biases (cf. Figure 13 and Table 4 in Shahmoradi & Nemiroff 2010).

Before the LGRB world model of Equation (1) is fit to BATSE observational data, it is necessary to convolve the model with the BATSE trigger threshold (Equation (5)). The study of the BATSE detection efficiency is well documented in a series of articles by the BATSE team (e.g., Pendleton et al. 1995, 1998; Paciesas et al. 1999; Hakkila et al. 2003; cf. Shahmoradi & Nemiroff 2011 for further discussion and references). Here, based on the observation that almost all 1366 BATSE LGRBs have durations of T₉₀ > 1 s, the primary trigger timescale for BATSE LGRBs is assumed to be 1024 ms. This eliminates the relatively complex dependence of the detection probability (η in Equation (5)) on the duration of the events. The probability of detection for an LGRB is then modeled by the cumulative density function (CDF) of log-normal distribution,

$$\begin{eqnarray} &&\eta ({\rm detection}|\mu _{{\rm thresh}},\sigma _{{\rm thresh}},L_{{\rm iso}},E_{p,z},z)\nonumber\\ &&\quad = \frac{1}{2} + \frac{1}{2}{\rm erf}\left (\frac{{\log {(P(L_{{\rm iso}},E_{p,z},z))} - \mu _{{\rm thresh}}}}{\sqrt{2}\sigma _{{\rm thresh}}} \right),\qquad \end{eqnarray} \tag{ B1 }$$

where P(L_iso, E_{p, z}, z) is the 1 s peak photon flux in the BATSE nominal detection energy range: 50–300 keV, and μ_thresh and σ_thresh are the detection threshold parameters to be determined by the model. The link between the 1 s peak photon flux (P) and the LGRB rest-frame variables (L_iso, E_{p, z}) and redshift (z) is provided by fitting a smoothly broken power law known as the Band model (Band et al. 1993) to the differential photon spectra of LGRBs,

$$\begin{equation} \hbox{{$\Phi (E) \propto \left\lbrace \begin{array}{@{}ll@{}}E^{\alpha } \exp (-{(1+z)(2+\alpha)E}/{E_{p,z}}) & E\le [E_{p,z}/(1+z)]\\ &\quad \times [(\alpha -\beta)/(2+\alpha)] \\ \ms E^{\beta } & \mbox{otherwise,} \end{array}\right.$}} \end{equation} \tag{ B2 }$$

such that

$$\begin{equation} P(L_{{\rm iso}},E_{p,z},z)=\frac{L_{{\rm iso}}}{4\pi {D_L}^2(z)} \frac{\int _{50}^{300} \Phi {\it dE}}{\int _{0.1/(1+z)}^{20000/(1+z)} E \Phi {\it dE}}. \end{equation} \tag{ B3 }$$

It has been shown by B10 that fixing the high-energy and low-energy photon indices of the Band model (Equation (B2)) to the corresponding population average α = −1.1, β = −2.3 produces only a negligible error of <0.05 dex in the resulting flux estimates. Given the uncertainties in the BATSE LGRB variables, in particular E_p estimates, such an assumption provides reasonable approximation for the calculation of the peak fluxes. A more accurate treatment, however, would be to include possible weak correlations that are observed between the Band model photon indices and the spectral peak energies of the bursts (cf. Shahmoradi & Nemiroff (2010, 2011) for a discussion of the correlations and potential origins).

The goodness of the log-normal CDF assumption for BATSE detection efficiency can be checked by a comparison of the resulting model predictions with the 1366 BATSE LGRBs' distribution of peak fluxes (Figure 2). The best-fit model prediction of BATSE trigger efficiency for the long-duration class of GRBs is compared to the nominal trigger efficiency of the BATSE 4B catalog for the class of soft-long bursts in Figure 10 (left panel). Although the difference between the two curves is significant, it does not necessarily imply a contradiction, given the fact that different methodologies and GRB samples were used to derive the two efficiency curves.³

To obtain the joint posterior for the unknown parameters of the LGRB world model of Equation (1) given BATSE data, the likelihood function of the model must be, in principle, constructed by correctly accounting for uncertainties in observational data (e.g., Eddington 1913; Jeffreys 1938). In addition, it is known that astronomical surveys at low S/Ns close to survey threshold can be potentially biased (e.g., Hogg & Turner 1998). A Bayesian multilevel methodology (e.g., Hobson et al. 2010) can incorporate the above corrections required to construct the likelihood function: under the assumption of normality for the uncertainties of LGRB variables in the BATSE catalog, each LGRB event—denoted by O_i, standing for the ith LGRB observation—has the likelihood $\mathcal {L}_i$ of having the true parameters $\boldsymbol O_i\equiv [P_{{\rm bol},i},S_{{\rm bol},i},E_{p,i},T_{90,i}]$ that is described by a 4D Gaussian pdf:

$$\begin{equation} \mathcal {L}_i(\boldsymbol O_i|\boldsymbol \mu _{i,0},\Sigma _{i,0})\sim \mathcal {N}(\boldsymbol O_i|\boldsymbol \mu _{i,0},\Sigma _{i,0}), \end{equation} \tag{ C1 }$$

where the location (μ_{i, 0}) and the scale (Σ_{i, 0}) parameters of the pdf are to be determined by the model and individual event data from the BATSE GRB catalog. Given $\mathcal {L}_i$ and the LGRB world model (Equation (1)) convolved with BATSE trigger efficiency (Equation (5)), the full Poisson likelihood function can be written as

$$\begin{eqnarray} &&{\mathcal {L}}({\rm Data}|{\rm Model\, Parameters})\,{=}\,\mathcal {A}^N\exp \left (\!{-}\mathcal {A}\!\int _{\boldsymbol O {\rm space}}\!\mathcal {R}_{{\rm obs}}(\boldsymbol O) \mathop {{d}}\nolimits ^{}\!\boldsymbol O\!\right) \nonumber \\ &&\quad\times \prod _{i=1}^{N} \int _{\boldsymbol O {\rm space}} \mathcal {R}_{{\rm cosmic}}(\boldsymbol \mu _{i,0},\Sigma _{i,0},z_i|{\rm Model\, Parameters})\mathcal {L}_i(\boldsymbol O_i|\nonumber\\ &&\quad\times\boldsymbol \mu _{i,0},\Sigma _{i,0}) \mathop {{d}}\nolimits ^{}\!\boldsymbol O, \end{eqnarray} \tag{ C2 }$$

where N = 1366, and $\mathcal {A}$ is a factor that properly normalizes the cosmic and the observed rates ($\mathcal {R}_{{\rm cosmic}}$ and $\mathcal {R}_{{\rm obs}}$). The term $\mathcal {R}_{{\rm cosmic}}$ in Equation (C2) acts as a prior for $\mathcal {L}_i$. In the absence of knowledge of the prior (i.e., $\mathcal {R}_{{\rm cosmic}}$, as in the case here), the empirical Bayes approach can provide an alternative solution, in which an ad hoc estimate of the model parameters based on the observed data (excluding uncertainties) serves as the prior for same data at the second level of analysis. Calculation of the normalization factor $\mathcal {A}$ involves an integration of the LGRB world model over the 5D space of LGRB variables and redshift, with a complex integration limit set by the BATSE trigger efficiency modeled in Equation (B1), as a function of L_iso, E_{p, z}, and z. In addition, since almost no redshift information is available for BATSE catalog of GRBs, the probability for each LGRB has to be marginalized over redshift, for which a range of z ∈ [0.1, ∞] is considered. These integrations make the maximization of the likelihood for all unknown parameters an extremely difficult task, perhaps challenging current computational technologies. Moreover, it has been known that GRB fluences and durations are likely underestimated close to detection threshold due to the so-called fluence–duration bias (e.g., Hakkila et al. 2000b, 2003). Such bias makes the overestimation correction of the fluence and duration as prescribed by Hogg & Turner (1998) unjustified before the fluence–duration bias effects are well quantified. The algorithms for calculating peak fluxes, however, appear to result in less biased measurements—even down to detector threshold—with negligible uncertainties (e.g., Stern et al. 2001). Given the above lines of reasoning and the computational limitations, the uncertainties on three BATSE LGRB variables, P_bol, S_bol, and T₉₀, are excluded from the likelihood function (Equation (C2)). The uncertainties on the E_p estimates of BATSE LGRBs are, however, significant compared to the three former variables and must be incorporated in the calculations of the likelihood. Nevertheless, it was realized after likelihood maximization that the exclusion of the E_p uncertainties—by fixing E_p values to the "bisector" estimates of Shahmoradi & Nemiroff (2010)—results in only negligible (≪1σ) changes in the model's best-fit parameters. In general, the use of the bisector line of OLS regression lines (e.g., Isobe et al. 1990) for estimation purposes is unfavored due to lack of a maximum likelihood interpretation. The special case of BATSE LGRBs here, however, turns out to be an exception.

In addition to cosmological time-dilation correction, it is common practice to make an energy correction to the temporal variables of Swift GRBs (e.g., Gehrels et al. 2006; B10), such as T₉₀, when translating the variable from the observer frame to the rest frame. For BATSE LGRBs, this energy correction is likely negligible, given the fact that the T₉₀ durations are calculated based on the total photon counts in the BATSE LAD energy range (e.g., Kouveliotou et al. 1993; Fishman et al. 1994), 20–2000 keV, which can be practically considered as bolometric. Nevertheless, it is expected that such an energy correction, if needed, would slightly relax the strong correlation of the rest-frame duration (T_{90, z}) with the total isotropic emission (E_iso) and the peak isotropic luminosity (L_iso).

The joint posterior distribution of the model parameters is obtained by maximizing the likelihood function of Equation (C2) convolved with a non-informative uniform prior on the location parameters and the standard choice of Jeffreys prior on the scale parameters (Jeffreys 1946). In order to efficiently sample from the 16D posterior density function, a variant of the MCMC methods known as AMH-MCMC is employed (e.g., Haario et al. 2001). The choice of an adaptive (versus classical) MH algorithm is very important, since the model parameters exhibit strong covariance (Table 3). To reduce the simulation runtime, all algorithms including AMH-MCMC are implemented in Fortran, which is by far the fastest and most efficient programming language for many intensive scientific calculations and number crunching (Loh 2010). In addition, the numerical integration in the definition of the luminosity distance (Equation (3))—encountered on the order of ≳ 10⁹ times during the full MCMC sampling—is greatly simplified by the analytical approximation method of Wickramasinghe & Ukwatta (2010). Due to the intrinsic sequential character of MCMC sampling methods, the parallelization of simulation algorithms (on either shared or distributed memory architecture machines) is impractical or at best inefficient for a single Markov chain. Nevertheless, to increase the MCMC sample size and more importantly, to ensure convergence to the global (versus local) extremum, the chain is initiated simultaneously at 10 random starting points in the parameter space on a 12-core desktop CPU. In general, convergence and good mixing occur within the first few thousands of iterations (burn-in period) given a suitable initial guess for the covariance matrix of the proposal distribution—here chosen to be multivariate Gaussian. The resulting mean and 1σ standard deviations of the model parameters are tabulated in Table 2.

The prediction power and consistency of the presented LGRB world model—based on BATSE data—can be easily checked against observational data from current and future gamma-ray experiments, in particular the Fermi satellite. All it takes is to construct a Monte Carlo universe of LGRBs, based on the best-fit parameters of the LGRB world model in Table 2, and compare the outcome with observational data. Although straightforward, the steps for such simulation and comparison are summarized below.

• 1.  
  A random redshift for each simulated LGRB is drawn from the redshift distribution of Equation (4) with parameters taken from Table 2. It is recommended to repeatedly randomly draw the set of model parameters from the full Markov chain samples⁴ instead of fixing the parameters to the mean values reported in Table 2.
• 2.  
  The four LGRB variables: L_iso, E_iso, E_{p, z}, T_{90, z} are randomly drawn from a 4D log-normal distribution with location (μ) and scale (i.e., the covariance matrix: Σ) parameters constructed from fitting results in Table 2. This can be easily and quickly done by noting that a multivariate log-normal distribution is equivalent to a multivariate Gaussian distribution in the logarithmic space of the above variables,

  $$\begin{equation} \log \boldsymbol O \equiv (\log (L_{{\rm iso}}),\log (E_{{\rm iso}}),\log (E_{p,z}),\log (T_{90,z})), \end{equation} \tag{ D1 }$$

  such that the 4D log-normal density function, $\mathcal {LN}$, of Equation (1) can be exactly replaced by a 4D Gaussian distribution,

  $$\begin{eqnarray} \mathcal {N}(\log \boldsymbol O|\boldsymbol \mu,\Sigma ) &=& (2\pi)^{-k/2}|\Sigma |^{{-1}/{2}}\nonumber\\ & \times & \exp (-(\log \boldsymbol O - \boldsymbol \mu)^{\prime }\Sigma ^{-1}(\log \boldsymbol O - \boldsymbol \mu)/2),\nonumber\\ \end{eqnarray} \tag{ D2 }$$

  for which the cosmic LGRB differential rate ($\mathcal {R}_{{\rm cosmic}}$) of Equation (1) will be

  $$\begin{equation} \mathcal {R}_{{\rm cosmic}} = \frac{{\it dN}}{d\log (L_{{\rm iso}}) d\log (E_{{\rm iso}}) d\log (E_{p,z}) d\log (T_{90,z}) dz}, \end{equation} \tag{ D3 }$$
• 3.  
  The above Monte Carlo universe of LGRBs can be then measured according to the instrument's detection efficiency (Section 2.2 and Equation (5)). For the case of BATSE LAD detectors, the trigger efficiency can be modeled according to the prescription in Appendix B.