SELECTION EFFECTS IN GAMMA-RAY BURST CORRELATIONS: CONSEQUENCES ON THE RATIO BETWEEN GAMMA-RAY BURST AND STAR FORMATION RATES

Gamma-ray bursts (GRBs) are the farthest sources, seen up to redshift z = 9.46 (Cucchiara et al. 2011), and if emitting isotropically they are also the most powerful (with E_iso ⩽ 10⁵⁴ erg s⁻¹) objects in the universe. Notwithstanding the variety of their peculiarities, some common features may be identified by looking at their light curves. A crucial breakthrough in this area has been the observation of GRBs by the Swift satellite, launched in 2004. With the instruments on board, the Burst Alert Telescope (15–150 keV), the X-Ray Telescope (0.3–10 keV), and the Ultra-Violet/Optical Telescope (170–650 nm), Swift provides a rapid follow-up of the afterglows in several wavelengths with better coverage than previous missions. Swift observations have revealed a more complex behavior of the light curve afterglow (O'Brien et al. 2006; Sakamoto et al. 2007) that can be divided into two, three, and even more segments in the afterglows. The second segment, when it is flat, is called plateau emission. A significant step forward in determining common features in the afterglow light curves was made by fitting them with an analytical expression (Willingale et al. 2007, hereafter W07). This provides the opportunity to look for universal features that could provide a redshift-independent measure of the distance from the GRB, as in studies of correlations between the GRB isotropic energy and peak photon energy of the νF_ν spectrum, E_iso − E_peak, (Lloyd & Petrosian 1999; Amati et al. 2009), the beamed total energy E_γ − E_peak (Ghirlanda et al. 2004; Ghirlanda et al. 2006), the luminosity-variability (L-V) (Norris et al. 2000; Fenimore & Ramirez-Ruiz 2000), L-E_peak (Yonetoku et al. 2004), and luminosity-τ lag (Schaefer 2003).

Dainotti et al. (2008, 2010), using the W07 phenomenological law for the light curves of long GRBs, discovered a formal anti-correlation between the X-ray luminosity at the end of the plateau L_X and the rest-frame plateau end-time, $T^*_a = T^{{\rm obs}}_a/ (1 + z)$, where $T^*_a$ is in seconds and L_X is in erg s⁻¹. The normalization and the slope parameters a and b are constants obtained by the D'Agostini fitting method (D'Agostini 2005). Dainotti et al. (2011a) attempted to use the L_X–$T^*_a$ (LT) correlation as a possible redshift estimator, but the paucity of the data and the scatter prevented them from a obtaining a definite conclusion, at least for a sample of 62 GRBs. In addition, a further step to better understand the role of the plateau emission has been made with the discovery of new significant correlations between L_X, and the mean luminosities of the prompt emission, <L_{γ, prompt} > (Dainotti et al. (2011b)). The LT anticorrelation is also a useful test for theoretical models such as the accretion models, (Cannizzo & Gehrels 2009; Cannizzo et al. 2011), the magnetar models (Dall'Osso et al. 2011; Bernardini et al. 2012; Bernardini 2012; Rowlinson et al. 2010; Rowlinson et al. 2013; Rowlinson et al. 2014), the prior emission model (Yamazaki 2009), the unified GRB and active galactic nucleus (AGN) model (Nemmen et al. 2012), and the fireshell model (Izzo et al. 2013). Moreover, Hascoet et al. (2014) and Van Eerten (2014b) consider both the LT and the L_X-<L_{γ, prompt} > correlation to discriminate among several models proposed for the origin of the plateau. In Leventis et al. (2014) and Van Eerten (2014a), a smooth energy injection through the reverse shock has been presented as a plausible explanation for the origin of the LT correlation. Furthermore, other authors were also able to reproduce and use the LT correlation to extend it in the optical band (Ghisellini et al. 2009), to extrapolate it into correlations of the prompt emission (Sultana et al. 2012), and to use the same methodology to build an analogous correlation in the prompt (Qi & Lu 2012). Finally, it has been applied as a cosmological tool (Cardone et al. 2009, 2010; Dainotti et al. 2013a; Postnikov et al. 2014). The impacts of detector thresholds on cosmological standard candles have also been considered (Shahmoradi & Nemiroff 2009; Petrosian et al. 1998; Petrosian 1999, 2002; Cabrera et al. 2007). However, because of large dispersion (Butler et al. 2010; Yu et al. 2009) and the absence of good calibration, none of these correlations allow the use of GRBs as good standard candles, as has been done, e.g., with Type Ia supernovae. An important statistical technique to study selection effects for treating data truncation in GRB correlations is the Efron & Petrosian (1992) method. Another way to study the same problem in GRB correlations, derived by modeling the high-energy properties of GRBs, has been reported in Butler et al. (2010). In the latter paper, it has been shown that well-known examples of these correlations have common features indicative of strong contamination by selection effects. We compare this procedure with the method introduced by Efron & Petrosian (1992) and applied it to the LT correlation (Dainotti et al. 2013b). The paper is organized as follows: Section 2 introduces the relation between the GRB rate (RGRB) and the star formation rate (SFR), Section 2.1 is dedicated to the analysis of a GRB scaling relation; in particular, we consider the LT correlation as an example, but the procedure described can be adopted for any other correlation. In Section 3 we describe how to build the GRB samples. In Section 4 we analyze the redshift evolution of the slope and normalization of the LT correlation. In Section 5 we study the selection effects related to simulated samples assuming different normalization and slope values. Then, in Section 6 we draw conclusions on the intrinsic slope of the LT correlation and on the evaluation of the redshift-dependent ratio between the RGRB and the SFR.

In order to understand the relation between GRBs and star formation, it is often assumed that the RGRB is proportional to the SFR, and the predicted distribution of the GRB redshift is compared to the observed distribution (Totani 1997; Mao & Mo 1998; Wijers et al. 1998; Porciani & Madau 2001; Natarajan et al. 2005; Jakobsson et al. 2006; Daigne & Mochkovitch 2007; Le & Dermer 2007; Coward 2007; Mao 2010). However, this relationship is not an easy task to handle because some studies show that GRBs do not seem to trace star formation in an unbiased manner (Lloyd & Petrosian 1999). Namely, the ratio between the RGRB and the SFR (RGRB/SFR) significantly increases with redshift (Kistler et al. 2008; Yuksel & Kistler 2007). This means that GRBs are more frequent for a given star formation rate density at earlier times. In fact, while observations consistently show that the comoving rate density of star formation is nearly constant in the interval 1 ⩽ z ⩽ 4 (Hopkins & Beacom 2006), the comoving rate density of GRBs appears to be evolving distinctly. In our approach we explicitly take into consideration this issue when we fit the observed RGRB with the model. Selection effects involved in a GRB sample are of two kinds: GRB detection and localization, and the redshift determination through spectroscopy and photometry of the GRB afterglow or the host galaxy. These problems have been objects of extensive study in the literature (Bloom 2003; Fiore et al. 2007; Guetta & Della Valle 2007). Moreover, the Swift trigger is very complex and the sensitivity of the detector is very difficult to parameterize exactly (Band 2006), but in this case, not dealing with prompt peak energy, we do not have to take into consideration the double truncation present in the data (Lloyd & Petrosian 1999). In the case of the plateau it is easier, since an effective luminosity threshold appears to be present in the data, which can be approximated by a 0.3–10 keV energy flux limit F_lim ≡ 2 × 10¹² erg cm² s¹ (Dainotti et al. 2013b). The luminosity threshold is then $L_{{\rm lim}} = 4\pi D_L^2(z,\Omega _M,H_0) F_{{\rm lim}}$, where D_L is the luminosity distance to the burst. Throughout the paper, we assume a flat universe with Ω_M = 0.28, Ω_λ = 0.72, and H₀ = 70 km s¹ Mpc¹. In our approach below several models are considered and then the one that best matches the RGRB with the SFR has been chosen.

2.1. GRB World Model

We derive a model capable of reproducing the observed Swift RGRB as a function of redshift, luminosity, and time of the plateau emission.

Rest-frame time and luminosity at the end of the GRB plateau emission show strong correlations as discovered by Dainotti et al. (2008) and later updated by Dainotti et al. (2010, 2011a, 2011b, 2013b). Therefore, all these quantities must be considered in deriving reliable rates. We characterize the RGRB as a product of terms involving the redshift z of the bursts, the isotropic equivalent luminosity release (0.3–10 keV) L_X, and the duration $T^{*}_a$.

Let us assume that a scaling relation exists so that the luminosity L_X(T_a) for a GRB with timescale $T^*_a$ at redshift z is given by

$$\begin{equation} \lambda =\alpha _0+\alpha _{\tau }\tau + \alpha _{\zeta }\zeta \end{equation} \tag{ 1 }$$

where we have introduced the compact notation

$$\begin{equation} \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\lambda = log \ L_X(T^*_a) \\ \ms \tau =log\ [T^*_a/(1+z)] \\ \ms \zeta =log\ (1+z). \end{array}\right. \end{equation} \tag{ 2 }$$

and the term ζ accounts for redshift evolution. The luminosity is normalized by the unit of 1 erg s⁻¹ and the time by the unit of 1 s, so that non-dimensional quantities are considered. All the observables in this model are computed in the rest frame because we are testing the role played by selection effects in the rest frame, the LT correlation rest frame being corrected. Independently, on the physical interpretation of this relation (in fact, there are several models that can reproduce it as we have mentioned in the Introduction), we can nevertheless expect GRBs to follow Equation (1) with a scatter of σ_λ. Moreover, the zero point α₀ may be known only up to a given uncertainty σ_α. Following the approach of Butler et al. (2010), applied for prompt correlations, we assume that λ can be approximated by a Gaussian distribution with mean λ₀, expressed in Equation (1), and the variance σ_int is the intrinsic scatter of the correlation. We also write the probability that a GRB with given (τ, ζ) values has a luminosity λ as follows:

$$\begin{equation} P_{\lambda }(\lambda,\tau,\zeta) \propto exp \left[-\frac{1}{2} \left[\frac{\lambda -(\alpha _0+\alpha _{\tau }\tau + \alpha _{\zeta }\zeta)}{\sigma _{\rm int}} \right]^2 \right] \end{equation} \tag{ 3 }$$

where $\sigma _{\rm int}^2 = \sigma _{\lambda }^2 + \sigma _{\alpha }^2$, σ_α is the uncertainty of the α₀ value, and σ_λ is the uncertainty on the luminosity value.

The approximation of a Gaussian distribution both for luminosity and time is motivated by the goodness of the fit, which gives a probability of P = 0.46 and P = 0.61, respectively, see Figure 1 and 2. We note that the mean (indicated with <>) <T_a > = 3.35 (s) with a variance $\sigma _{T_a}=0.77$ (s) and <L_X > = 48.04 (erg s⁻¹) with a variance $\sigma _{L_X}=1.37$ (erg s⁻¹) are represented, respectively, in Figures 1 and 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to obtain the number of GRBs with a given luminosity λ, we need to integrate over the distributions of τ and ζ. We will assume, for simplicity, that τ follows a truncated Gaussian law.

$$\begin{equation} P_{\tau }(\tau) \propto \left\lbrace \begin{array}{@{}l@{\quad }l@{}}exp[-\frac{1}{2}(\frac{\tau -\tau _0}{\sigma _{\tau }})^2] & \tau _L<\tau <\tau _U \\ 0 & \tau \le \tau _L \ or \ \tau \ge \tau _U \end{array}\right. \end{equation} \tag{ 4 }$$

where τ_L and τ_U respectively indicate the lower limit and upper limit of the observed τ distribution and τ₀ is the mean value of this distribution. The limits of τ are taken from an updated sample of $T^*_a$ composed of 176 GRB afterglows, with firm redshift determination, from 2005 January until 2014 July. The analysis follows the criteria adopted in Dainotti et al. (2013b).

If we assume that the GRBs trace the cosmic SFR, we can model their redshift distribution following Butler et al. (2010):

$$\begin{equation} P_z(z) \propto \frac{\dot{\rho }_*(z)}{1+z}\frac{dV}{dz} \end{equation} \tag{ 5 }$$

where $\dot{\rho }_{*}(z)$ is the comoving RGRB density, V is the universal volume, and the factor (1 + z) accounts for the cosmic time dilatation and

$$\begin{equation} \frac{dV}{dz} \propto \frac{r^2(z)}{E(z)} \end{equation} \tag{ 6 }$$

with r(z) being the comoving distance and E(z) = H(z)/H₀ being the Hubble parameter normalized to its present-day value.

Collecting the different terms, we can finally write the true, detector-independent event $\mathcal {N}$ differential rate, for z, $\log T^*_a$, and log L_X:

$$\begin{equation} \frac{d\mathcal {N}}{d\lambda d\tau dz} \propto \Psi (z) P_{\lambda }(\lambda, \tau, \zeta) P_{\tau }(\tau)P_z(z). \end{equation} \tag{ 7 }$$

We note here that we have introduced the term of the evolution in redshift, Ψ(z) = (1 + z)^α, following the approach of Lloyd & Petrosian (1999); Dermer (2007), and Robertson & Ellis (2012). In Dermer (2007), assuming that the emission properties of GRBs do not change with time, they find that the Swift data can only be fitted if the comoving rate density of GRB sources exhibits positive evolution to z > 3–5. In our approach, we introduce evolution starting from z ⩾ 0.99.

Therefore, using the above expression for P_τ, we find that the number of GRBs with luminosity in the range (λ, λ + dλ) and redshift between z and z+dz is

$$\begin{eqnarray} \frac{d\mathcal {N}}{d\lambda dz} &\propto & \Psi (z) \frac{\dot{\rho }_*(z)(dV/dz)}{1+z} \frac{\mathcal {F}_{\tau _U}-\mathcal {F}_{\tau _L}}{\sqrt{8\pi \sigma ^2_{\tau }}} \nonumber\\ && \times exp \left[-\frac{1}{2} \left[\frac{ \lambda -(\alpha _0+\alpha _{\tau }\tau + \alpha _{\zeta }\zeta)}{\sqrt{\sigma ^2_{\rm int}}} \right]^2 \right] \end{eqnarray} \tag{ 8 }$$

where $\mathcal {F}_{\tau _U}$ and $\mathcal {F}_{\tau _L}$ are the respective error functions⁷ of the lower and upper limits of the time distribution. Note that Equation (8) is defined up to an overall normalization constant, which can be solved by imposing that the integral of $d\mathcal {N}/d\lambda dz$ over (λ, z) gives the total number of observed GRBs. Actually, this is not known since we do not observe all GRBs but only those passing a given set of selection criteria. However, we will only be interested in the fraction of GRBs in a cell in the two-dimensional (2D) (λ, z) space so that we do not need this quantity.

We are aware that we do not map out the true LT relation, given selection effects and the observed LT relation. Doing this would require modeling the selection of the GRB sample itself (using the gamma-ray threshold) and also seeking to understand the tie between the GRB flux and the afterglow L_X. However, the relation between the flux and L_X has already been studied by Dainotti et al. (2013b) and was reported briefly in the previous section above. Here we compute the new limit related to the updated sample, as shown in the middle panel of Figure 3.

Zoom In Zoom Out Reset image size

The RGRB given by Equation (8) has been derived by implicitly assuming that all the GRBs can be detected notwithstanding their observable properties. This is actually not the case. As an example, we will consider hereafter the LT correlation, although the formalism and the method we will develop can be easily extended to whatever scaling law. For the LT case, there are two possible selection effects. First, each detector has an efficiency that is not the same for all the luminosities. Only GRBs with λ > λ_L, where λ_L is the lowest detectable luminosity for a given instrument, can be detected while all the GRBs with λ larger than a threshold luminosity λ_U will be found.

Moreover, it is likely that the efficiency of the detector is not constant, but is rather a function of the luminosity. We will therefore introduce an efficiency function $\mathcal {E}_{\lambda }(\lambda)$ whose functional expression is not known in advance, but can only take values in the range of (0,1). A second selection effect is related to the time duration of the GRB. Indeed, in order to be included in the sample used to calibrate the LT correlation, the GRB afterglow has to be measured over a sufficiently long timescale to make it possible to fit the data and extract the relevant quantities. If τ is too small (as shown in Dainotti et al. 2013b; the minimum rest-frame time is 14 s), few points will be available for the fit, while, on the contrary, large τ values will give rise to afterglow light curves, which could be well sampled by the data. Again, we can parameterize these effects by introducing a second efficiency function $\mathcal {E}_{\tau }(\tau)$ so that the final observable rate is the following:

$$\begin{equation} \frac{d\mathcal {N}_{\rm obs}}{d\lambda dz} \propto \frac{d\mathcal {N}}{d\lambda dz}\times \mathcal {E}_{\lambda }(\lambda)\mathcal {E}_{\tau }(\tau). \end{equation} \tag{ 10 }$$

We point out that our formulation, which takes into account the efficiency functions $\mathcal {E}_{\lambda }(\lambda)$ and $\mathcal {E}_{\lambda }(\tau)$ in the final observed RGRB, is similar to the approach by Robertson & Ellis (2012) in Equation (1), in which the additional factor K is presented. K is equivalent to our $\mathcal {E}_{\lambda }(\lambda)$ and $\mathcal {E}_{\lambda }(\tau)$.

It is worth noting that Equation (10) is actually still a simplified description. Indeed, it is in principle possible that other selection effects take place involving observable quantities not considered here, as for example β and the redshift. However, these parameters enter into the determination of λ so that one can (at least in a first-order approximation) convert selection cuts on them in a single efficiency function, depending only on λ (for the dependence of the flux on the redshift, see the left panel of Figure 3). However, as we can see from Figure 3 β is constant with redshift, and there is no correlation between those two quantities; in fact, the Spearman correlation coefficient is ρ = −0.062. Nevertheless, Equation (10) provides a reasonably accurate description of the observable RGRB.

In order to evaluate Equation (10), there are different quantities to determine. First, we need to set the scaling coefficients (α₀, α_τ, α_ζ) and the intrinsic scatter σ_int. Second, the mean and variance of the τ distribution (τ₀, σ_τ) has to be given. Finally, an expression for the cosmic SFR $\dot{\rho }_*(z)$ has to be assigned. None of these quantities is actually available. In principle, one could assume an SFR law and fit for the model parameters to a large enough GRB sample with measured (λ, τ, ζ) values. To this end, one should know the selection function $\mathcal {E}_{\lambda }(\lambda)\mathcal {E}_{\tau }(\tau)$, which is not the case. Studies of how light curves would appear to a gamma-ray detector here on Earth have been performed (Kocevski & Petrosian 2013). In this paper the prompt emission pulses are investigated and the conclusion is that even a perfect detector that observes over a limited energy range would not faithfully measure the expected time dilation effects on a GRB pulse as a function of redshift.

Nevertheless, here we study detector threshold effects on afterglow properties. Our aim is to investigate how the ignorance of the efficiency function biases the estimate of the correlation coefficients. We can therefore rely on simulated samples based on a realistic intrinsic rate. We proceed as schematically outlined below.

• 1.  
  We assume that the available data represent reasonably well the intrinsic τ distribution so that we can infer (τ₀, σ_τ) from the data themselves. We set τ_{L, U} = τ₀ ± 5σ_τ, thus symmetrically cutting the Gaussian distribution at its extreme ends.
• 2.  
  Based on the shape of the cosmic SFR (Hopkins & Beacom 2006), we assume a broken power law for the comoving RGRB density:

  $$\begin{equation} \dot{\rho }_*(z) \propto \left\lbrace \begin{array}{@{}l@{\quad }l@{}}(1+z)^{g_0} & z\le z_0 \\ \ms (1+z)^{g_1} & z_0\le z \le z_1 \\ \ms (1+z)^{g_2} & z \ge z_1 \end{array}\right. \end{equation} \tag{ 11 }$$

  where the relative normalizations are set so that $\dot{\rho }_*(z)$ is continuous at z₀ = 0.97 and z₁ and (z₀, z₁) = (0.97, 4.50), (g₀, g₁, g₂) = (3.4, −0.3, −8.0). Moreover, besides Equation (11), we employed other shapes of the SFR (Li 2008, Robertson & Ellis 2012; Kistler et al. 2013) to obtain the observed RGRB density. The one used by (Li 2008) is

  $$\begin{equation} \dot{\rho }(z)=a + b\times Log(1 + z). \end{equation} \tag{ 12 }$$

  The a and b parameters are

  $$\begin{equation} (a, b) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}(-1.70, 3.30) & z \le 0.993 \\ \ms (-0.727, 0.0549) & 0.993 \le z \le 3.80 \\ \ms (2.35,-4.46) & z \ge 3.80 \\ \end{array}\right. \end{equation} \tag{ 13 }$$

  Robertson & Ellis (2012) defined the SFR as

  $$\begin{equation} \dot{\rho }(z) = \frac{a+b(z/c)^f}{1+(z/c)^d}+g, \end{equation} \tag{ 14 }$$

  where they have a = 0.009 M_☉ yr⁻¹ Mpc⁻³, b = 0.27 M_☉ yr⁻¹ Mpc⁻³, c = 3.7, d = 7.4, and g = 10⁻³ M_☉ yr⁻¹ Mpc⁻³.Instead, Kistler et al. (2013) defined the SFR as

  $$\begin{equation} \dot{\rho }(z) = \dot{\rho }_0 \times \left[(1 + z)^{a \psi } + \left(\frac{1 + z}{B} \right)^{b \psi } + \left(\frac{1 + z}{C} \right)^{c \psi } \right]^{\frac{1}{\psi }}, \end{equation} \tag{ 15 }$$

  with slopes of a = 3.4, b = −0.3, and c = −2.5, and breaks at z₁ = 1 and z₂ = 4 corresponding to $B = (1+z_1)^{1-\frac{a}{b}} \sim 5160$ and $C=(1+z_1)^{\frac{(b-a)}{c}}\times (1+z_2)^{\frac{(1-b)}{c}} \sim 11.5$, respectively, and ψ = −10.Finally, we compare the fitted functions obtained with these four methods with our data distribution. The most reliable fit for our parameters is the SFR used by Li (2008); see Figure 4 where the best fit among linear (upper panel) and polynomial (lower panel) (λ) functions are considered. Moreover, we adopted the constraints for the redshift-dependent ratio between SFR and RGRB adopted by Robertson & Ellis (2012). In this paper a modest evolution (e.g.,Ψ(z) ≈ (1 + z)^α) with −0.2 ⩽ α ⩽ 1.5, where the peak probability occurs for α ≈ 0.5 is consistent with the long GRB prompt data (P ≈ 0.9). These values can be explained if GRBs occur primarily in low-metallicity galaxies that are proportionally more numerous at earlier times. We note that in our approach we assumed no evolution at low redshift for z ⩽ 0.99 consistently with the posterior probability in Robertson & Ellis (2012) in which no evolution is possible at the 2-σ level. However, because a constant Ψ(z) is also ruled out (Robertson & Ellis 2012), we fit the normalization parameters and the evolution factors, obtaining Ψ(z) ≈ (1 + z)^−0.2 for 0.993 ⩽ z ⩽ 3.8 and Ψ(z) ≈ (1 + z)^0.5 for z ⩾ 3.8. These values of the evolution are compatible with Robertson & Ellis (2012). Regarding the observed RGRB, we obtained that the best efficiency functions are possible both for two polynomial and two linear functions, as we show in Figure 4. Tables 1 and 2 show the probability that the density rate matches the afterglow plateau RGRB, assuming those efficiency functions.
• 3.  
  For given (α₀, α_τ, α_ζ, σ_int) values, we divide the 2D space (λ, z) in $\mathcal {M}$ cells and, for each cell, compute the fraction of GRBs in it:

  $$\begin{equation} f_{\rm sim}(\lambda _i,z_i)=\frac{\int _{\lambda _i-\Delta \lambda }^{\lambda _i+\Delta \lambda }d\lambda \int _{z_i-\Delta z}^{z_i+\Delta z}dz \frac{d\mathcal {N}}{d\lambda dz}}{\int _{\lambda _{\rm min}}^{\lambda _{\rm max}}d\lambda \int _{z_{\rm min}}^{z_{\rm max}}dz \frac{d\mathcal {N}}{d\lambda dz}} \end{equation} \tag{ 16 }$$

  where we set

  $$\begin{equation} (\lambda _{\rm min}, \lambda _{\rm max}) = (42.0, 52.0),(z_{\rm min}, z_{\rm max})=(0,10). \end{equation} \tag{ 17 }$$

  We find it more efficient to change the variable from z to ζ when dividing the 2D space in 10 × 10 square cells.
• 4.  
  For each given cell, we generate $\mathcal {N}_{ij}=f_{\rm sim}(\lambda _i, \zeta _j) \times \mathcal {N}_{\rm sim}$ GRBs (with $\mathcal {N}_{\rm sim}$ the total number of objects to simulate) by randomly sampling (λ, ζ) within the cell boundaries and computing τ by solving Equation (1).
• 5.  
  To take into account the selection effects, for each GRB we generate two random numbers (u_τ, u_λ), uniformly sampling the range (0, 1), and only retain the GRB if $u_{\tau }\le \mathcal {E}_{\tau }(\tau)$ and $u_{\lambda }\le \mathcal {E}_{\lambda }(\lambda)$. Note that, as a consequence of this cut, the final number $\mathcal {N}_{\rm obs}$ of observed GRBs is smaller than the input one $\mathcal {N}_{\rm sim}$.
• 6.  
  Finally, for each one of the $\mathcal {N}_{\rm obs}$ selected GRBs, we generate new (τ_obs, λ_obs) values, extracting from Gaussian distributions centered on the simulated (τ, λ) values and with a 1% variance. We also assign to each GRB an error set in such a way as to be similar to what is actually obtained for GRBs with comparable (τ, λ) values.

Zoom In Zoom Out Reset image size

Table 1. Efficiency Function Parameters for the Power-law $\mathcal {E}_{\lambda }$ and no cut on τ, i.e., $\mathcal {E}_{\tau }=1$.

Id	λL	λU	$\mathcal {E}_{\lambda }$	PGRB, rate
PL1	44.34	50.86	1.25	⩽10−4
PL2	43.64	49.87	2.99	0.003
PL3	43.77	50.74	1.65	0.53
PL4	44.77	49.59	2.04	0.001
PL5	44.14	50.83	0.23	0.54

Notes. P_{GRB, rate} is the goodness of fit between our data and the observed RGRB density and thus how well the data fit the observed GRB density rate. To compute the probability, we compute the χ² test, which performs a statistical hypothesis test in which the sampling distribution of the test statistic is a χ² distribution when the null hypothesis is true, in order to determine whether or not there is a significant difference between the expected frequencies and the observed frequencies.

Download table as:  ASCIITypeset image

Table 2. Same as Table 1 but for the Polynomial Functions

Id	λL	λU	$\mathcal {E}_1$	$\mathcal {E}_2$	$\mathcal {E}_3$	$\mathcal {E}_4$	PGRB, rate
PoL1	44.90	49.14	0.46	0.01	0.24	0.80	0.54
PoL2	41.10	50.23	0.60	0.95	0.05	0.53	⩽10−4
PoL3	43.57	49.09	0.71	0.79	0.07	0.34	0.019
PoL4	44.37	49.52	0.51	0.03	0.46	0.78	0.15
PoL5	43.03	50.06	0.79	0.36	0.63	0.40	0.001

Note. P_{GRB, rate} is the goodness of fit between our data and the observed RGRB density.

Download table as:  ASCIITypeset image

The above procedure allows us to build a simulated GRB sample taking into account both the intrinsic properties of any scaling relation and the selection effects induced by the instrumental setup. Moreover, we have referred to an actual GRB sample in order to set both the limits on τ, ζ, and λ, and the typical measurement errors. Therefore, we can rely on these simulated samples to investigate the impact of selection effects on the recovered slope and intrinsic scatter of the given correlation. To this end, the last ingredient we need is a functional expression for the efficiency functions. Since these are largely unknown, we are forced to make some arbitrary guess. Therefore, we consider two different cases. First, we assume that there is no selection on τ, i.e., we set $\mathcal {E}_{\tau } = 1$. Two functional expressions are then used for $\mathcal {E}_{\lambda }$, namely a power law:

$$\begin{equation} \mathcal {E}_{\lambda }(\lambda) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}0 & \lambda <\lambda _L \\ (\frac{\lambda -\lambda _L}{\lambda _U-\lambda _L})^{\mathcal {E}_{\lambda }} & \lambda _L\le \lambda \le \lambda _U \\ 1 & \lambda >\lambda _U \end{array}\right. \end{equation} \tag{ 18 }$$

and a fourth-order polynomial, i.e.,

$$\begin{equation} \mathcal {E}_{\lambda }(\lambda) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}0 & \lambda <\lambda _L \\ \frac{\mathcal {E}_1\tilde{\lambda }+\mathcal {E}_2\tilde{\lambda }^2+\mathcal {E}_3\tilde{\lambda }^3+\mathcal {E}_4\tilde{\lambda }^4}{\mathcal {E}_1+\mathcal {E}_2+\mathcal {E}_3+\mathcal {E}_4} & \lambda _L\le \lambda \le \lambda _U \\ 1 & \lambda >\lambda _U \end{array}\right. \end{equation} \tag{ 19 }$$

with $\tilde{\lambda }=(\lambda -\lambda _L)/(\lambda _U-\lambda _L)$. We try different arbitrary choices for the parameters entering both expressions of $\mathcal {E}_{\lambda }$ in order to investigate to what extent the results depend on the exact choice of the efficiency function; see Figure 5. In a second step, we abandon the assumption $\mathcal {E}_{\tau } = 1$, to assume for it the same functional expression used for $\mathcal {E}_{\lambda }$, with the same choices for the parameters, but with different upper and lower limits depending on τ_U and τ_L; see Figure 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As we have already mentioned in the previous paragraph the polynomial and the linear model for the (λ) are unknown, and therefore assumptions need to be made. We chose these forms because both the normalization and the slope of the LT correlation depend on the redshift either with a polynomial or with a simple power law. Therefore, these choices for the selection functions take into account this redshift dependence. Namely, we consider a model redshift dependence because of the corresponding dependence of both luminosity and time. This has been already shown in Dainotti et al. (2013b) and is also shown in the middle panel of Figure 3 for the updated data sample. To study the behavior of the redshift evolution we plot the slope and the normalization values versus the redshift. These are obtained from the average values for the data set divided into 5 bins (see Figure 7) and into 12 bins (see Figure 8). As we can see from both Figures 7 and 8 that the normalization parameter α₀ decreases as the redshift increases, while the slope parameter α_τ shows the opposite trend. Goodness of fit is given by the probability P = 0.79 for the data set divided into 5 bins and P = 0.87 for the one divided into 12 bins for the linear case, while for the polynomial model it is P = 0.99 and P = 0.94 for the data set divided in 5 and 12 bins, respectively. These results show that both polynomial and linear fits are possible.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The simulated samples generated as described above are input into the same Bayesian fitting procedure we use with real data. For each input (α_τ, α_ζ, α₀, σ_int) parameter, we simulate a ∼50 GRB sample setting $\mathcal {N}_{\rm sim} = 200$, while the number of observed GRBs depends on the efficiency function used. We fit these samples assuming no redshift evolution in Equation (1), i.e., forcing α_ζ = 0 into the fit so that, for each simulated sample, the fitting procedure returns both the best fit and the median and 68% confidence range of the parameters (α_τ, α₀, σ_int). In order to investigate whether or not the selection effects impact the recovery of the input scaling laws, we fit linear relations of the form:

$$\begin{equation} x_f = ax_{\rm inp} + b \end{equation} \tag{ 20 }$$

where x_inp is the input value and x_f can be either the best fit (denoted as x_bf) or the median x_fit value. When fitting the above linear relation, we use the χ² minimization for x_bf, while a weighted fit is performed for x_fit with weights $\omega _i = 1/\sigma _i^2$ where σ_i is the symmetrized 1σ error. Note that the label i here runs over the simulations performed for each given efficiency function.

5.1. No Redshift Evolution

Here, we consider input models with α_ζ = 0, i.e., no redshift evolution of the scaling law (1). It is worth noting that such an assumption is actually well motivated since it has been demonstrated in Dainotti et al. (2013b) that luminosity is almost unaffected by redshift evolution, while time undergoes redshift evolution for high redshift only. From our point of view, however, this case allows us to directly quantify the impact of the efficiency functions on the recovery of the scaling correlation parameters, since any deviation will only be due to the selection effects and not due to any attempt of compensating the missed evolution with z.

5.1.1. No Selection on τ ($\mathcal {E}_{\tau } = {\it 1}$)

We start by considering the idealized case of no selection of τ, i.e., we force $\mathcal {E}_{\tau } = 1$, and set the $\mathcal {E}_{\lambda }$ parameters as listed in Tables 1 and 2 for the power-law and polynomial expressions, respectively. As an example, Figure 9 shows the results for the efficiency function, while Table 3 summarizes the (a, b) coefficients of the linear fit between the input and recovered quantities. The closer a is to 1, the less the parameter is biased, while b ≠ 0 should not be taken as evidence for bias. This result is in perfect agreement with the intrinsic correlation slope, which is $-1.07^{+0.09}_{-0.14}$ (Dainotti et al. 2013b), when we consider the selection functions that return values of a closer to 1 as the best choices. If Equation (20) is fulfilled, we can estimate the relative bias:

$$\begin{equation} \frac{\Delta x}{x}=\frac{x_{\rm inp}-x_f}{x_{\rm inp}}=1-a-\frac{b}{x_{\rm inp}}, \end{equation} \tag{ 21 }$$

Zoom In Zoom Out Reset image size

Table 3. Slope a and Zero Point b of the Fitted versus Input Parameters for Both the Best-fit and Median Values (Labeled with Subscripts bf and fit, Respectively)

Id	$(a,b)^{\tau }_{bf}$	$(a,b)^{\tau }_{\rm fit}$	$(a,b)^{\alpha _0}_{bf}$	$(a,b)^{\alpha _0}_{\rm fit}$	$(a,b)^{\sigma }_{bf}$	$(a,b)^{\sigma }_{\rm fit}$	Δx/x
PL1	(0.953, 0.010)	(0.959, 0.013)	(0.928, 3.688)	(1.000, −0.073)	(0.593, 0.354)	(0.616, 0.355)	0.004
PL2	(0.914, −0.008)	(0.873, −0.024)	(1.013, −0.836)	(0.989, 0.292)	(0.689, 0.299)	(0.643, 0.341)	0.002
PL3	(0.880, −0.052)	(0.937, −0.013)	(0.965, 1.729)	(1.008, −0.513)	(0.683, 0.340)	(0.664, 0.352)	0.003
PL4	(0.946, 0.024)	(0.964, 0.024)	(0.995, −0.076)	(1.086, −4.905)	(0.614, 0.364)	(0.585, 0.380)	0.006
PL5	(0.916, −0.030)	(0.962, 0.004)	(1.033, −2.067)	(0.828, 9.095)	(0.716, 0.333)	(0.679, 0.356)	0.005

Notes. The upper script denotes the parameter fitted with τ, α₀, and σ referring to α_τ, α₀, and σ_int, respectively. Δ_x/x is the bias for each efficiency function considered.

Download table as:  ASCIITypeset image

so that we can accept b ≠ 0 if x_inp is much larger than b. This is indeed the case for x_inp = α₀, which takes typical values (∼50) that are much larger than the b ones in Table 3.

From the proximity between solid and dashed lines, which respectively represent the best-fit line and the no-bias line when x_inp = x_f, in the corresponding panels of Figure 9, we see that, for the power-law efficiency function (and no cut on τ), both the slope and the zero point of the scaling relation are correctly recovered. The reason why is that the relative bias is negligibly small, notwithstanding the values of the parameters setting $\mathcal {E}_{\lambda }$. This is particularly true if one relies on the median values as estimates since they are typically consistent with the no-bias line within less than 2σ.

The above results have been obtained considering a power-law $\mathcal {E}_{\lambda }$ so that it is worth investigating whether or not they critically depend on this assumption. We have therefore repeated the analysis for the polynomial $\mathcal {E}_{\lambda }$ models in Table 2, obtaining the results in Table 4. A comparison with the values in Table 4 shows that the (a, b) coefficients are similar so that one could preliminarily conclude that the shape of the efficiency function does not play a major role in the determination of the bias. Actually, although the functional expressions are different, both the power-law and the polynomial selection functions are qualitatively similar with $\mathcal {E}_{\lambda }$ increasing with λ over a comparable range. Although such a behavior is likely common to any reasonable $\mathcal {E}_{\lambda }$, we can not exclude a priori that non-monotonic selection functions do actually exist. What the results would be in such a case is not clear so that we prefer to be cautious and conclude that the bias is roughly the same whichever monotonic $\mathcal {E}_{\lambda }(\lambda)$ function is used, but not for all the possible $\mathcal {E}_{\lambda }$ functions. For the non-monotonic shape of the selection function, see Stern et al. (2001), in which an assumed detection efficiency function, defined as the ratio of the number of detected test bursts to the number of test bursts applied to the data versus the expected peak count rate, is given by

$$\begin{equation} E(c_e)=0.70 \times \left[1-exp \left[-\left(\frac{c_e}{c_{e,0}} \right)^2 \right] \right]^{\nu }, \end{equation} \tag{ 22 }$$

Table 4. Same as Table 3 but for the Polynomial $\mathcal {E}_{\lambda }$ Model

Id	$(a,b)^{\tau }_{bf}$	$(a,b)^{\tau }_{\rm fit}$	$(a,b)^{\alpha _0}_{bf}$	$(a,b)^{\alpha _0}_{\rm fit}$	$(a,b)^{\sigma }_{bf}$	$(a,b)^{\sigma }_{\rm fit}$	Δx/x
PoL1	(0.950, 0.020)	(0.928, 0.011)	(1.099, −5.545)	(1.178, −9.935)	(0.647, 0.350)	(0.673, 0.349)	0.004
PoL2	(1.095, 0.128)	(1.075, 0.105)	(1.030, −1.662)	(1.023, −1.278)	(0.681, 0.337)	(0.625, 0.372)	0.0008
PoL3	(0.984, 0.063)	(0.936, 0.015)	(0.711, 15.258)	(0.988, 0.376)	(0.741, 0.289)	(0.685, 0.336)	0.007
PoL4	(0.870, −0.025)	(0.969, 0.052)	(0.730, 14.103)	(0.734, 13.872)	(0.630, 0.351)	(0.582, 0.381)	0.009
PoL5	(1.004, 0.069)	(0.972, 0.030)	(0.963, 11.772)	(1.002, −0.374)	(0.581, 0.384)	(0.549, 0.402)	0.18

Download table as:  ASCIITypeset image

where c_{e, 0} = 0.097 counts s⁻¹ cm⁻² and ν = 2.34 are two constants. However, quoting from Stern et al. (2001), the best possible efficiency quality has still not yet been achieved because in fact the detection efficiency depends on the peak count rate rather then on the time-integrated signal.

5.1.2. Selection Cuts on Both τ and λ

We now consider the case where the total selection function may be factorized as $\mathcal {E}(\tau, \lambda) = \mathcal {E}_{\tau }(\tau)\mathcal {E}_{\lambda }(\lambda)$ with both $\mathcal {E}_x(x)$ functions being given by power-law or fourth-order polynomial expressions. We consider 10 different arbitrary choices for both cases. Note that we have to increase $\mathcal {N}_{\rm sim}$ to 300 in order to have $\mathcal {N}_{\rm obs} = 80\hbox{--}100$ as for the models discussed in the previous subsection.

Table 5 gives the (a,b) coefficients for the different models considered. A comparison with Table 3 shows that, on average, the bias on the parameters is roughly the same with the median values giving smaller deviations and significant bias on σ_int only. A more detailed analysis, however, shows that, in the $\mathcal {E}_{\tau }=1$ case, biases larger than 5% are of the order of 10%. Namely, from Table 3 and 4 we show that the relative biases, Δ_x/x, both in the linear and the polynomial case, give very small values from 0.2% to 0.9%, with the only exception of 1 polynomial function, in which the bias is 18%, thus giving a P = 10% of having larger bias than 5%. If we consider selection cuts on both τ and λ we notice in Table 5 that the number of the bias whose value is greater than 5% is four, thus increasing their probability of occurring (P ∼ 40%). This can be qualitatively explained by noting that a cut on λ only removes the points in the luminosity axis, thus possibly shifting the best-fit relation, but not changing the slope. On the contrary, removing points also along the horizontal τ axis can change the slope α_τ too and hence also affects (α₀, σ_int) because of the correlation among these parameters and α_τ. Similar results are obtained when both $\mathcal {E}_x$ functions are modeled with fourth-order polynomials, so we will not discuss this case here. We stress that when Δ_x/x is larger than 6%, the slope of the correlation is farther from −1.0 compared to cases in which Δ_x/x ⩽ 0.06. In fact, in the first case the slope values range from 0.86 to 0.91; see the functions PLTa4 and PLTa6 in Table 5. These values are not compatible in 1 σ with the claimed intrinsic slope of the correlation, $-1.07^{+0.09}_{-0.14}$. If we consider, instead, the lowest Δ_x/x, then we obtain ranges of a = (0.94, 0.99), thus showing full compatibility in 1 σ with the intrinsic slope. In this way we have quantitatively confirmed the existence of the $L_X-T^{*}_a$ correlation with the same intrinsic slope as in Dainotti et al. (2013b) if appropriate selection functions are chosen.

Table 5. Same as Table 3 but for the Selection Cuts on both (τ, λ) and Power-law $\mathcal {E}_x$ Functions

Id	$(a,b)^{\tau }_{bf}$	$(a,b)^{\tau }_{\rm fit}$	$(a,b)^{\alpha _0}_{bf}$	$(a,b)^{\alpha _0}_{\rm fit}$	$(a,b)^{\sigma }_{bf}$	$(a,b)^{\sigma }_{\rm fit}$	Δx/x
PLTa1	(0.842, −0.080)	(0.920, −0.016)	(1.099, −5.615)	(1.064, −3.650)	(0.615, 0.363)	(0.585, 0.387)	0.04
PLTa2	(0.972, −0.016)	(0.963, −0.027)	(1.112, −5.980)	(1.148, −7.746)	(0.720, 0.306)	(0.674, 0.335)	0.005
PLTa3	(0.960, −0.025)	(0.974, −0.005)	(0.985, 0.791)	(1.022, −1.228)	(0.725, 0.329)	(0.690, 0.345)	0.004
PLTa4	(0.877, −0.021)	(0.933, −0.005)	(1.005, −0.611)	(0.959, 2.196)	(0.702, 0.304)	(0.644, 0.345)	0.09
PLTa5	(0.860, −0.056)	(0.901, −0.026)	(0.911, 4.592)	(0.838, 8.655)	(0.687, 0.338)	(0.613, 0.382)	0.06
PLTa6	(0.904, −0.010)	(0.919, −0.003)	(0.734, 14.060)	(0.397, 32.119)	(0.678, 0.326)	(0.680, 0.330)	0.08
PLTa7	(0.915, −0.044)	(0.998, 0.035)	(1.119, −6.577)	(1.119, −6.474)	(0.736, 0.311)	(0.722, 0.327)	0.02
PLTa8	(0.967, −0.003)	(0.961, −0.017)	(1.027, −1.485)	(1.178, −9.515)	(0.731, 0.287)	(0.705, 0.313)	0.03
PLTa9	(0.933, −0.005)	(0.904, −0.031)	(1.026, −1.530)	(0.804, 10.289)	(0.641, 0.357)	(0.670, 0.352)	0.06
PLTa10	(1.040, 0.056)	(0.997, 0.028)	(0.753, 13.308)	(0.822, 9.345)	(0.736, 0.329)	(0.713, 0.344)	0.04

Download table as:  ASCIITypeset image

Here we built a general method to evaluate selection effects for GRB correlations, not knowing a priori the efficiency function of the detector used. We have tested this method on the LT correlation. We chose a set of GRBs and assumed Gaussian distributions for the variables involved, for luminosity and time, and also a particular shape for the RGRB density. We simulated a mock sample of data in order to consider the selection effects of the detectors. As we can see in Section 3, assuming the correct observed RGRB density shape was not an easy task. In fact, we explored different methods (Li 2008; Robertson & Ellis 2012; Kistler et al. 2013; Hopkins & Beacom 2006) that use several SFR shapes to understand which one best matches the afterglow plateau data distribution, including the selection functions; see Figure 4. The most reliable fit for the GRB plateau data is the SFR used by Li (2008), while the best efficiency functions for (λ) that match the GRB density rate can be both two polynomial and two linear functions; see Figure 4. Tables 1 and 2 show the probability that the density rate fits the afterglow plateau RGRB assuming those efficiency functions. However, we assumed that there could be selection effects for both luminosity and time. In particular, the bias is roughly the same whichever monotonic efficiency function for the luminosity detection $\mathcal {E}_{\lambda }$ is taken. In Tables 3 and 4 we show that the relative biases, Δ_x/x, both in the linear and the polynomial cases, give very small values from 0.2% to 0.9%, with the only exception being the 1 polynomial function, in which the bias is 18%, thus giving P = 10% of having a larger bias than 5%. If we consider selection cuts on both τ and λ, we notice in Table 5 that the number with bias whose value is greater than 5% is four, thus increasing the probability of having such biases (P ∼ 40%). In addition, we studied selection effects in the LT correlation assuming also a combination of the luminosity and time detection efficiency functions. Different values for the parameters of the efficiency functions in the detectors are taken into account as described in Section 5. This gives distinct fit values that are inserted into Equation (20) and allow us to study the scattering of the correlation and its selection effects. We have quantitatively confirmed the existence of the $L_X-T^{*}_a$ correlation with the same intrinsic slope as in Dainotti et al. (2013b) if appropriate selection functions are chosen. In particular, when Δ_x/x is larger than 6%, the slope of the correlation is farther from −1.0 compared to cases in which $\frac{\Delta _x}{x} \le 0.06$. The lowest Δ_x/x leads to ranges of a = (0.94, 0.99), thus showing full compatibility in 1 σ with the intrinsic slope. Finally, the fact that the correlation is not generated by the biases themselves is a significant and further step toward considering a set of GRBs as standard candles and its possible and useful application as a cosmological tool.

This work made use of data supplied by the UK Swift Science Data Centre at the University of Leicester. We are particularly grateful to V. F. Cardone for the initial contribution to this work. M.G.D. and S.N. are grateful to the iTHES Group discussions at Riken. M.D. is grateful to the support from the JSPS Foundation (No. 25.03786). N.S. is grateful to JSPS (No.24.02022, No.25.03018, No.25610056, No.26287056) & MEXT(No.26105521), and R.D.V. is grateful to 2012/04/A/ST9/00083.