Accounting for Selection Bias and Redshift Evolution in GRB Radio Afterglow Data

Gamma-ray Bursts (GRBs) are highly energetic events that can be observed at extremely high redshift. However, inherent bias in GRB data due to selection effects and redshift evolution can significantly skew any subsequent analysis. We correct for important variables related to the GRB emission, such as the burst duration, $T_{90}^*$, the prompt isotropic energy, $E_{\rm iso}$, the rest-frame end time of the plateau emission, $T_{\rm a, radio}^*$, and its correspondent luminosity $L_{\rm a, radio}$, for radio afterglow. In particular, we use the Efron-Petrosian method presented in 1992 for the correction of our variables of interest. Specifically, we correct $E_{\rm iso}$ and $T_{90}^*$ for 80 GRBs, and $L_{\rm a, radio}$ and $T_{\rm a, radio}^*$ for a subsample of 18 GRBs that present a plateau-like flattening in their light curve. Upon application of this method, we find strong evolution with redshift in most variables, particularly in $L_{\rm a, radio}$, with values similar to those found in past and current literature in radio, X-ray and optical wavelengths, indicating that these variables are susceptible to observational bias. This analysis emphasizes the necessity of correcting observational data for evolutionary effects to obtain the intrinsic behavior of correlations to use them as discriminators among the most plausible theoretical models and as reliable cosmological tools.


Introduction
Gamma-ray bursts (GRBs) can be observed at extremely high redshift, making them a valuable tool for studies of the early Universe. The ability to observe these highly energetic objects at redshifts out to z = 9.4 [1][2][3][4][5][6][7][8][9][10][11] has created significant interest in using them as standardizable candles, similar to Type Ia supernovae. However, observations of GRBs have shown a very diverse population with few common characteristics.
Phenomenologically, GRBs are characterized by the main event, called the prompt emission, which is usually observed in gamma-rays, hard X-rays and sometimes in optical, while the afterglow is the counterpart in soft X-rays (≈66% of observed GRBs), in optical (≈38% of observed GRBs) and sometimes in radio (≈6.6% of observed GRBs). GRB radio afterglows are very difficult to observe, indeed, similar to the X-ray observations which are characterized by the detector limits, and additional difficulties rise due to the limited allocated time for the follow-up observations in the radio band after the GRB trigger. Bursts are classified following the duration of the prompt episode (T 90 ). The population of short GRBs (sGRBs) usually has harder spectra and a duration of less than 2 s. In contrast, the population of long GRBs (lGRBs) has softer spectra and a duration larger than 2 s [12]. However, this classification is still in debate, and so in some bursts, it is not clear if GRBs with intermediate features belong to short or long GRBs [13,14].
LGRBs are associated with the core collapses of dying massive stars [15,16] and sGRBs by merging two compact objects; a black hole (BH) with a neutron star (NS) and two NSs [17][18][19].
A crucial breakthrough in the analysis of GRB features is the discovery of the plateau emission by the Neils Gehrels Swift Observatory [20]. The plateau emission is a flat part of the lightcurves which follows the decay phase of the prompt emission [21][22][23][24]. In the current paper we focus on properties resulting from this plateau emission, as well as both the prompt and afterglow emission. In general, attempts have been made to find standardizable properties, such as a plateau of GRBs or through prompt and afterglow correlation studies. We here mention a few of them: Yonetoku et al. [25], Dai and Wang [26], Ghirlanda [27], Dainotti et al. [28], Amati et al. [29], Dainotti et al. [30][31][32][33][34], Dainotti [35].
However, it is clear from past and current studies [36][37][38][39] that observations of GRBs are further susceptible to selection bias and evolutionary effects, which may change the results of any subsequent analysis and can substantially impact the results related to cosmological application of GRB relations [40]. In GRB studies, it is, therefore, crucial to know whether the studied correlations are intrinsic or artificially created as a result of observational biases and redshift evolution. "Redshift evolution" is the dependence of the variable of interest on redshift, and thus "independent of redshift" indicates the absence of such evolution.
In the study of GRB correlations, all variables must be computed in the rest-frame, as we are comparing objects at different epochs. This introduces another source of redshift dependence included in the definition of luminosity distance: where H 0 is the Hubble constant at the present day and Ω M is the matter density in a flat Universe assuming the equation of state parameter w = −1. Indeed, usually one of the variables in the correlation is either a luminosity or energy which, by definition, depends on the luminosity distance. Ideally, all correlations we use must be corrected for redshift evolution, if any, requiring the removal of any existing redshift dependence. There do exist statistical techniques that are capable of correcting for these effects, as well as correcting for data truncation from detector limits [41][42][43]. Among the methods to remove evolution, we consider here the Efron-Petrosian (EP) [41] method. The EP method is a well-established example of these kinds of techniques, and has been used to recover intrinsic relationships in many correlations in the past [7,40,[44][45][46][47][48].
Lloyd et al. [44] discuss the correlation between E p , or the peak of the νF ν spectrum, with flux and fluence in GRBs, later investigated in the rest-frame and known as the E p − E iso relation [49]. A further modification of this relation is the one discovered by Yonetoku et al. [25] in which E p is correlated with the prompt isotropic luminosity, L iso . The EP method provides an explanation on how to perform analysis on truncated data, and in Yonetoku et al. [25], Amati et al. [49] it is illustrated that the method is capable of recovering the correlation present in the original "parent" sample with the truncated data.
This technique has been further explored regarding the luminosity function and formation rate of sGRBs in a recent study by Dainotti et al. [48]. They look at the intrinsic distributions of these variables using the EP method, and introduce a method of accounting for incompleteness of redshift data with the Kolmogorov-Smirnov (KS) test (this is described in more detail in Section 2.2). They find a strong evolution of luminosity with redshift, em-phasizing the necessity of this correction. The analysis presented in Dainotti et al. [48] is also relevant, as it emphasizes that both sGRBs and lGRBs undergo strong redshift evolution.
It should also be noted that though this method is mainly applied to GRB correlation studies, it has been also successfully applied in studies of Active Galactic Nuclei as well [45].
Among GRB correlations in particular we focus our attention to the rest-frame time at the end of the plateau emission, T * a,radio , and its correspondent luminosity L a,radio , this is an extension in radio of the so-called 2D Dainotti relation in X-rays [29][30][31] and optical [50]. For the very recent analysis on the 2D Dainotti relation in radio see Levine et al. (2021) in preparation. For a review of the subject of GRB correlations and selection biases in the prompt and afterglow see Dai and Wang [26], Ghirlanda [27], Dainotti [35], Dainotti and Del Vecchio [51], Dainotti et al. [52], Dainotti and Amati [53].
One of the main problems in the application of GRB relationships as theoretical model discriminators and as cosmological tools is the fact that correlations must be intrinsic to the physics and not induced by biases. There are several examples of how the correlations are used to interpret theoretical models both in the prompt and afterglow emission. The photospheric emission and the Comptonization models [54][55][56][57][58][59][60][61][62] are the two main models used to test the E peak -E iso and the Yonetoku et al. [25] correlations, the latter between E peak and the isotropic energy in the prompt emission. Otherwise, the parameter space pinpointed by those correlations can be the effect of selection biases and not of the true underlying physics. To this end, it is necessary to apply these correlations as model discriminators only after correction for such biases. Indeed, for example the plateau emission in X-rays and optical, which reconciles with the existence of the 2D Dainotti relation, can be derived through the equations of a fast rotating NS, the so-called magnetar model [17][18][19][63][64][65]. In Rowlinson et al. [66], Rea et al. [67], Stratta et al. [68] the derivation of the parameter space of the magnetic field and spin period have been computed accounting for selection biases and redshift evolution. The current status in the literature is that only a few correlations have been corrected for selection biases and evolutionary effects through the EP method, such as Dainotti et al. [7,48,69,70,71].
Specifically, Dainotti et al. [7] examine this correlation in X-ray for a sample of 101 GRBs that present a plateau, or flattening, in their light curves. After correction for evolutionary effects using the EP method, they conclude that the observed correlation is intrinsic at the 12 σ level. In mimicking the evolution of each variable with redshift, they tested both a simple and more complex model, finding similar results in both cases. Dainotti et al. [40] further examine the importance of these corrections when studying the cosmological properties of GRBs, applying the EP method to a simulated correlation between luminosity and time at the end of the plateau emission for 101 GRBs and testing whether a 5 σ deviation from the intrinsic values strongly changes the cosmological results. They demonstrated that their results change with this deviation by 13% regarding the values of Ω M , emphasizing the necessity of applying such corrections. The problem of evolution of the variables and their correction is not only important for GRB-cosmology studies, but also for more general cosmological studies. Indeed, in Dainotti et al. [72] it has been shown that there is indication of a possible evolution even of the Hubble constant. If this is not due to selection biases, a new physical cosmological model which relies on alternative theories, such as the modified theory of gravity, must be accounted for.
Regarding the correlations in GRB afterglows, Lloyd-Ronning et al. [46,73] discussed the correlation of not only luminosity with redshift, but also isotropic energy, E iso , T * 90 , and the jet opening angle, θ j , for a sample of 376 GRBs. They emphasize the difficulty of obtaining intrinsic values for these quantities due to inherent biases in observation methods, and additional truncation from detector limits that can introduce false correlations in the data [74]. They find strong evolution with redshift for each of these variables, indicating that achromatic properties of GRBs are also susceptible to selection bias. A further study by Lloyd-Ronning et al. [47] discusses the evolution of θ j with redshift in greater detail, using the EP method to recover the intrinsic behavior of the jet opening angle.
In this study, we seek to determine whether the strong evolution of E iso and T * 90 vs. redshift initially found by Lloyd-Ronning et al. [46,73] is still the same for GRBs with observed radio afterglow. In addition to the isotropic energy, we apply the EP method to the luminosity, and break time in radio wavelengths to determine if these variables are strongly affected by inherent bias and evolutionary effects.
We here point out that we are aware that the plateau sample is a subsample of a more extended population of plateaus that we cannot see. We have fixed the issue of the biases related to the redshift evolution and due to the selection threshold with the Efron-Petrosian method; however, we cannot account for the missing population of GRBs for which the followup has not been tackled. Nonetheless, it is crucial to discuss the L a,radio versus the redshift, since this correlation has been studied in X-rays and optical extensively and it is important to investigate if this correlation holds true in radio with comparable slopes to optical and X-ray. The first step to investigate the correlation is to determine if the variables involved are subjected to redshift evolution and selection biases.
In summary, the main point of the paper is the study of the redshift evolution and the removal of selection biases through the EP method. The analysis of the true correlations can be done only if we first determine the evolution among the variables. The plateau emission has been extensively investigated in X-rays and optical, but so far there has not been a statistical analysis of the existence of the plateau in radio. The radio observations of the plateau emission can cast a light on whether or not the end time of the plateau is indeed a jet break. This point can be revealed only through such a study. The evaluation of the jet break allows one to better understand the evolution of the GRB and its physics in relation to the standard fireball model [75,76] or other scenarios. This paper is organized as follows: in Section 2, we discuss the selection of our sample, as well as the formulation of the EP method and its application to our sample. In Section 3, we present the results of this analysis, and in Section 4, we discuss the implications of our study, as well as a comparison to previous studies, and present our conclusions.

Variables of Interest
The EP method uses a modified version of Kendall's τ statistics to test for independence of variables in a truncated data set. Here, τ is defined as: with R i defined as the rank E i = (1/2)(i + 1) defined as the expectation value, and V i = (1/12)(i 2 + 1) defined as the variance. For a more complete discussion of the method and the algebra involved, see Dainotti et al. [7,34], Efron and Petrosian [41], Singal et al. [45], Dainotti et al. [77], Petrosian et al. [78], Lloyd-Ronning and Petrosian [79]. Here, we use the EP method to determine the impact of redshift evolution and selection bias on four variables: T * 90 , where the star denotes the rest-frame, E iso , the radio light curve break time T * a,radio , and the radio luminosity at the time of break L a,radio . These variables are considered to be they are pertinent to the correlations analyzed in Levine et al. (2021 in preparation). Throughout our analysis, we consider these variables in logarithmic scale for convenience.
We look at the log E iso and log T * 90 for a sample of 80 GRBs with observed radio afterglow published in the literature [80][81][82][83][84][85][86][87][88][89][90][91][92][93][94][95][96][97][98]. Values of log E iso are taken from the literature. If no E iso value could be found, the log E iso , in units of ergs, is calculated using the equation: where S is the fluence in units of erg cm −2 , D 2 L (z) is the luminosity distance assuming a flat ΛCDM model with Ω M = 0.3 and H 0 = 70 km s −1 Mpc −1 (see Equation (1)), and K is the correction for cosmic expansion [99]: with β as the spectral index of the GRB. Fluence and β values are taken from the literature.
To analyze the impact of selection bias and redshift evolution for log L a,radio and log T * a,radio , we first fit each of the 80 GRBs with a broken power law (BPL) according to the formulation: where F a refers to the flux at the break, T * a refers to the rest-frame time of break, α 1 refers to the slope before the break, and α 2 refers to the slope after the break. We can only obtain values of log L a,radio and log T * a,radio for light curves that show a "plateau" feature, or a flattening of the light curve before a clear break. In our analysis, we consider a light curve to display a plateau if |α 1 | < 0.5. Therefore, we discard fits to light curves with scattered observations, unreliable error bars, or shapes incompatible with a BPL and plateau. In our subsequent analysis we include those light curves whose ∆χ 2 analysis of the BPL best-fit parameters are suitable following the Avni [100] methodology. After the rejection process, we are left with 18 GRBs that present a plateau and clear break in the light curve.
The luminosity log L a,radio in units of erg s −1 is computed at time log T * a,radio using the equation: where F a is the observed flux at T a,radio , D 2 L (z) is defined as in Equation (3), and K is the k-correction: with β as the radio spectral index and α 1 as the fitted BPL temporal index before the break. β values are taken from Chandra and Frail [80] or other literature-if no β value could be found, the average of published spectral indices, β = 0.902 ± 0.17, was assigned. We show the distribution of spectral indices for the plateau sample in Figure 1.

Limiting Fluxes, Fluences, and Times
The EP method can overcome selection bias for a particular variable of interest, but it must first be determined if the variable is dependent or independent of redshift.
It is then necessary to define limiting values for each of the variables. In Dainotti et al. [7], it was demonstrated that a good choice of limiting times and luminosities retains at least 90% of the total sample. For time variables, log T * 90 and log T * a,radio (in units of seconds), we define a general form for the limiting values as T min (1+z) , where T min is the minimum observed time. We need to choose a compromise between a limit which is representative of the population of data points, but still retains most of the sample size. It has been shown in Monte Carlo simulations in [7] that such a strategy with limiting values is accurate. For log T * 90 , we find the best limiting duration to be log T * 90 min,obs = −0.54 s, with a limiting boundary defined as −0.30 (1+z) s, which excludes 5/80 (<10%) GRBs. The limiting line for log T * 90 is shifted at a higher value to be allow the sample data to be representative of the whole population. For log T * a,radio , we find the limiting time to be the observed minimum log T * a,radio min,obs s, which does not exclude any data points. For the isotropic energy, we instead define the limiting energy according to the methodology of Dainotti et al. [7], in which the limiting fluence should be representative of the population while including at least 90% of the sample. We use the following formula: where S lim is the fluence limit. For our sample, we define S lim as 6.3 × 10 −8 erg cm −2 .
Applying this limit excludes 8/80 GRBs, which is 10% of our sample. In all the method described here we use GRBs that have log E iso > log E iso,lim , log T * 90 > log T * 90 , log T a,radio > log T a,radio,lim , and log L a,radio > log L a,radio,lim . For the luminosity, however, a caveat should be posed when we consider the total distribution of the parent population of GRBs with and without redshift (see [48]).
Using the method presented in Dainotti et al. [48], we compare the parent sample of all GRBs with observed radio afterglow and known peak flux to a smaller "subsample" of GRBs with known peak flux and known redshift. We then apply cuts to the data by defining limiting fluxes at regular intervals. Considering only the data with values above the limiting fluxes, fluences and time, we conduct a two-sample Kolmogorov-Smirnov (KS) test between the data of the total sample and the data for which the limiting cuts have been applied to determine the distribution of the probability that the subsample was drawn from the parent sample, as well as the geometric distance between the two samples as determined by the KS test. We take the limiting flux to be the value of f lim where the probability as a function of limiting flux reaches a plateau in which the probability that two samples are drawn by the same population is 100%. In our sample, we find this limit to be log f lim = −17.2. We define the flux throughout our analysis in units of erg cm −2 s.
We show the distribution of the parent sample and subsample in the left panel of Figure  2, with the limiting line shown in red. We plot this probability as a function of flux limit (blue), as well as the distance between the distributions (orange), in the right panel of Figure 2.

Removing Selection Bias and Redshift Evolution
After defining the limiting lines, we can then mimic the evolution of this variable with redshift using a simpler function of redshift, f (z) = 1 (1+z) δ . We here adopt the choice of a simple function, but in principle it is possible to use a more complex function as already shown in [77] and obtain compatible results. Using a modified version of Kendall's τ statistics, we can compute an evolutionary function f (z) where the slope of the function is defined by δ. The values for δ where τ = 0 represent the removal of evolutionary effects. We define errors on δ out to 1 σ, corresponding to the δ values where τ = ±1.

Results
We here clarify that the purpose of our analysis is to show how similar our results are, compared to other ones in the literature given that our sample size differs from other studies for E iso and T * 90 (our sample has 80 GRBs), while this is the first time in the literature that we compare the results for the L a,radio and T * a,radio (our sample has 18 GRBs) with the previous results in the literature performed in X-rays and optical. This is an essential comparison to allow the determination of the intrinsic nature of the L a,radio and T * a,radio correlation and to check for the universality of the results related to the evolutionary functions for these variables with the EP method.
Using the procedure outlined above, we correct log E iso , log T * 90 using the formulation , where all quantities that have symbol are the variables for which the evolution has been removed, thus they are no longer dependent on the redshift. We find δ T * 90 = −0.65 ± 0.27, with the 1 σ errors defined as the average of the values of τ = 1 and τ = −1, and δ E iso = 0.39 ± 0.88. Figure 3 shows these results-the top left panel shows log T * 90 for the sample of 80 GRBs as a function of redshift, with the limiting value shown in red. The top right panel highlights the evolutionary function for log T * 90 , with dashed lines at τ = 0, ±1. The same plots for log E iso are shown in the bottom panels. For the plateau sample of 18 GRBs, we use the same formulation to obtain log L a,radio = log L a,radio − log((1 + z) δ L a,radio ) and log T * a,radio = log T * a,radio − log((1 + z) δ T * a,radio ). We find the values of δ for log L a,radio and log T * a,radio as δ T * a,radio = −1.94 ± 0.86 and δ L a,radio = 3.15 ± 1.65. These results are shown in Figure 4-log T * a,radio as a function of redshift is shown in the top left panel, with the limiting values in red. The evolutionary function is shown in the top right panel, with dashed blue lines at τ = 0, ±1. The bottom two panels display the same plots for log L a,radio .

Discussion and Conclusions
Using a sample of 80 GRBs with observed radio afterglow, we test the use of the EP method on log E iso and log T * 90 for a subsample 80 GRBs, and on log L a,radio and log T * a,radio for a subsample of 18 GRBs. We find that when considering log E iso and log T * 90 , we obtain indices for the parameter of redshift evolution, δ, of δ T * 90 = −0.65 ± 0.27 and δ E iso = 0.39 ± 0.88, while the values of δ for log L a,radio and log T * a,radio are log T * a,radio as δ T * a,radio = −1.94 ± 0.86 and δ L a,radio = 3.15 ± 1.65.
For log T * 90 , for log T a,radio , and log L a,radio , we find relatively strong evolution of each variable with redshift. The luminosity presents the strongest correlation with redshift, emphasizing the necessity of correction for these effects before using data in correlation analysis. E iso , by contrast, appears to be the most independent from redshift, with the smallest value for |δ|.
We find that our values are comparable to values obtained in previous studies. A study by Lloyd-Ronning et al. [46] of the cosmological evolution of isotropic energy E iso , burst duration T * 90 , jet opening angle θ j , and luminosity L j reports δ values compatible with our findings T * 90 and L a,radio within 1 σ, and values compatible with our δ E iso within approximately 2 σ. Specifically, they find a value of δ E iso = 2.3 ± 0.5, which agrees with our value of δ E iso = 0.39 ± 0.88 within 2.17 σ. We also see agreement with δ T * 90 = −0.8 ± 0.3, with a 0.55 σ difference from our value of δ T * 90 = −0.65 ± 0.27, and in the luminosity with δ L j = 3.5 ± 0.5, a 0.2 σ difference from our value of δ L a,radio = 3.15 ± 1.65.
These results also agree with previous values of δ for L a and T * a in X-ray and optical wavelengths. Dainotti et al. [7] conduct a similar analysis of the luminosity and rest-frame end time of the plateau emission using X-ray data. Their value for correction for δ T * a , reported as δ T * a,X = −0.85 ± 0.3, is compatible with our value of δ T * a,radio = −1.94 ± 0.86 within 1.23 σ.
However, they find a very slow evolution in luminosity, with a value of δ L a,X = −0.05 ± 0.35, which is a 1.94 σ difference from our value of δ L a,radio = 3.15 ± 1.65. This discrepancy is likely due in part to the small sample size, which may exaggerate the extent of the evolution present in our sample, but may also be due to differences in the behavior of the X-ray and radio emission.
We have corrected the luminosity and time in X-rays with 222 GRB lightcurves with a given redshift, presenting plateaus according the Willingale et al. [23] model and in optical with 181 GRBs with plateaus taken from Dainotti et al. [50], but with the additional analysis of 80 GRBs found in the literature. For tackling this analysis, we followed the same procedure described in the current paper. The results of this analysis reports σ, and δ L a,X = 2.42 ± 0.58, which is a 0.44 σ difference from our value. They also report values in optical of δ T * a,opt = −2.11 ± 0.49 and δ L a,opt = 3.96 ± 0.58, which both agree with our result within 1 σ.
In general, it can be seen that a larger sample size is preferred when applying the EP method. In our case, for the sample pertinent to E iso , and T * 90 we choose the limiting values while excluding ≤ 10% of the overall sample. However, for the smaller sample of 18 plateau GRBs, the limiting values are chosen so that we do not exclude any data points due to the small sample size. In addition, this conservative choice would allow us to have smaller error bars on the slope of the evolutionary functions. However, the δ values obtained for L a,radio and T * a,radio are similar to values found in previous studies of larger sample sizes, thus indicating that the EP method is still successful even with a small dataset.
GRB correlations in radio afterglows related to the plateau emission are crucial to understand if the jet break is coincident with the end of the plateau emission. To investigate this point, a multiwavelength analysis not only in optical and X-rays must be performed together with the radio data. Since the evolution of the variables with redshift can change the time at which the break happens, it is crucially important to correct for the redshift evolution. This study is also the preliminary step to the investigation of the intrinsic nature of the plateau emission correlations.
After our analysis on the radio observations both for the prompt emission in relation to the variables of E iso and T * 90 and for the afterglow emission in relation to L a,radio , T * a,radio , we can conclude:

1.
After testing intrinsic properties of a GRB, such as E iso and T * 90 , as well as properties such as L a,radio and T * a,radio , we see T * 90 , L a,radio and T * a,radio present strong correlation with redshift, thus indicating that they are susceptible to redshift evolution.

2.
The δ values obtained in this work agree with those of previous studies, indicating that this trend of strong correlation with redshift still holds true in radio wavelengths.

3.
The study of these evolutionary functions is the first necessary step to determine the true intrinsic nature of the correlations, object of our study.    [97] and Rhodes et al. [98].