Diagnosing Circumburst Environment with Multiband Gamma-Ray Burst Radio Afterglows

It has been widely recognized that gamma-ray burst (GRB) afterglows arise from interactions between GRB outflow and circumburst medium, while their evolution follows the behaviors of relativistic shock waves. Assuming the distribution of circumburst medium follows a general power-law form, that is, $n = A_{\ast} R^{-k}$, where $R$ denotes the distance from the burst, it is obvious that the value of density-distribution index $k$ can affect the behaviors of the afterglow. In this paper, we analyze the temporal and spectral behaviors of GRB radio afterglows with arbitrary $k$-values. In the radio band, a standard GRB afterglow produced by forward shock exhibits a late-time flux peak, and the relative peak fluxes as well as peak times at different frequencies show dependencies on $k$. Thus with multi-band radio peak observations, one can determine the density profile of circumburst medium by comparing the relations between peak flux/time and frequency at each observing band. Also, the effects of trans-relativistic shock waves, as well as jets in afterglows are discussed. By analyzing 31 long and 1 short GRBs with multi-band data of radio afterglows, we find that nearly half of them can be explained with uniform interstellar medium ($k=0$), $\sim 1/5$ can be constrained to exhibiting stellar wind environment ($k=2$), while less than $\sim 1/3$ samples show $0


INTRODUCTION
Gamma-ray bursts (GRBs) are stellar explosions arising from the collapses of massive stars (for long bursts, e.g., see Woosley 1993;Paczyński 1998;MacFadyen & Woosley 1999) or mergers of binary compact stars (for short bursts, e.g., see Paczynski 1986;Narayan et al. 1992;Gehrels et al. 2009;Abbott et al. 2017). The afterglows of GRBs are the results of interactions between relativistic ejecta from central engines and the circum-burst medium. The shock waves produced during such interactions can accelerate the swept-up electrons, thus giving rise to nonthermal radiation, such as synchrotron emission (e.g., see Piran 1999;van Paradijs et al. 2000;Mészáros et al. 2002). Our understanding of GRB afterglows has been greatly improved with the discoveries of GRB multiband afterglows (Costa et al. 1997;van Paradijs et al. 1997;Frail et al. 1997;Zhang 2007).
Assuming the profile of the circumburst density n follows a power-law form, that is, n = A * R −k , where R denotes the distance from GRB central engine, it is obvious that different values of the power-law index k lead to different light-curve behaviors (e.g., see Wu et al. 2005). The most widely discussed circumburst density model is the homogeneous interstellar medium (ISM) case with k = 0 (e.g., see Sari et al. 1998 & Kumar 2004), while the stellar-wind environment with k = 2 has also been proposed (e.g., see Dai & Lu 1998;Mészáros et al. 1998;Chevalier & Li 2000;Panaitescu & Kumar 2000, 2004Kobayashi & Zhang 2003;Wu et al. 2003Wu et al. , 2004Zou et al. 2005). Previously, it was thought that because the long GRBs originate from massive stars with significant mass losses via stellar winds before their demise, such events should preferably occur in stellar-wind environments. For example, Starling et al. (2008) performed a joint fitting for X-ray to IR afterglows of 10 GRBs, and put constraints on k for half of the sample, with 4 windlike cases and 1 ISM one. However, Panaitescu & Kumar (2002) found that half of their 10 GRBs favor a k = 0 medium; while Curran et al. (2009) made k constraints for 6 out of 10 GRBs, with 2 consistent with the ISM environment, 2 with the wind case, and another 2 samples compatible with both k = 0 and 2 cases. Yi et al. (2013) and van Eerten (2014) presented various relationships between observables in GRB X-ray and optical afterglows for arbitrary k. However, Yi et al. (2013) have drawn a typical value of k ∼ 1 for 19 Swift long bursts by analyzing their early afterglows that originated from forward-reverseshock interactions, implying a mass-loss scheme of the GRB progenitors other than stellar wind.
In order to give a better understanding of circumburst medium distributions, it is necessary to investigate this issue with methods other than early X-ray/Optical observations. In this work, we utilize a new way to constrain the density distribution of the circumburst medium, that is, to put a constraint on k by comparing multiband GRB radio afterglow peak times/fluxes with the corresponding frequencies. Observationally, the light curves of GRB radio afterglows usually show a double-peaked structure. In the rest frame, the typical peak times are 0.1 − 0.2 days and 2 days after triggers of prompt emissions, respectively (Chandra & Frail 2012). Generally speaking, the first peak can be explained with forward-reverse-shock interactions at early times, while the second peak is due to late-time forward-shock evolution and is unique for the radio band. In fact, because a strong dependency exists between GRB afterglow behaviors in the radio band and the underlying synchrotron radiation, Barniol Duran (2014) and Beniamini & van der Horst (2017) have already put constraints on the GRB microphysical parameters with GRB radio peaks, including the fraction of electron and magnetic energy among shock energy, that is, ε e and ε b , respectively.
In our circumburst density analysis, we ignore the effects of early radio peaks due to reverse shocks and only take the second peak into consideration. Since the peak flux F peak,ν and the peak time t peak,ν of this late peak depend on k, the value of k can be deduced with multiband radio-peak observations. Although in theory, the k value can also change the evolution of radio afterglows, it should be noted that radio observations are largely constrained by instrumental sensitivities, refractive interstellar scintillation (Chandra & Frail 2012), or the intrinsic properties of GRBs (Hancock et al. 2013).
Thus, the precise shapes of the complete light curves are usually hard to obtain. Hence, the main advantage of our radio-peak method is that the peak times/fluxes can be determined more easily and reliably due to the high fluxes of such peaks, as long as the afterglow light curves are sampled frequently enough. And because radio afterglows occur at a larger radius than early forwardreverse shocks, combining density profiles for earlier afterglows using the method proposed by Yi et al. (2013), and those for later radio peaks as shown in this work, a more complete circumburst medium profile at various R can be obtained.
This paper is organized as follows. In Section 2, we search for the peak time t peak,ν in the analytical light curve of the GRB radio afterglow with an arbitrary k, and calculate the relations between t peak,ν and F peak,ν with observing frequencies, with the role played by the k value clarified. Meanwhile, the effects of transrelativistic shocks, as well as jets, are discussed. In Section 3, we apply our theoretical predictions to 31 long GRBs, as well as 1 short burst, with multiband radio afterglow observations, and draw a conclusion that most of them can be explained with transrelativistic shock waves under ISM environments, which is different from Yi et al. (2013). This implies that the circumburst density profile may change with radius. In Section 4 our results are discussed and summarized. In our calculations, the ΛCDM cosmology is adopted, with cosmological parameters H 0 = 71km s −1 Mpc −1 , Ω m = 0.3, and Ω Λ = 0.7.

GRB RADIO LIGHT CURVES AND MULTIFREQUENCY PEAK EVOLUTION
In this section, our goal is to derive the relations between GRB radio afterglow peak times t peak,ν versus k, as well as peak fluxes F peak,ν versus k, i.e., the kdependent expressions of the power-law indices a and b for t peak,ν ∝ ν −a , and F peak,ν ∝ ν b , respectively. We investigate the analytical behaviors of multiband afterglow with arbitrary circumburst density-distribution index k in order to get the expressions of a and b. Ignoring the early afterglow evolution due to forward-reverseshock interactions, we consider the shock with radius R > R dec only, in which R dec is the deceleration radius of relativistic shock, with a corresponding deceleration time t dec = R dec (1 + z) /2Γ 2 0 c, a cosmological redshift of z, and an initial bulk Lorentz factor of Γ 0 . By this time, the GRB afterglow has already entered the selfsimilar phase. Here we do calculations for semiradiative shock waves for generality, although by the time of the emergence of radio peaks, i.e., several days after the initial GRB trigger, the shock waves are usually adiabatic. Assuming the circumburst density can be described as n = A * R −k = A (R/R * ) −k (k = 0 for uniformly distributed interstellar medium, k = 2 for stellar wind, with a typical value of the characteristic length R * ∼ 10 17 cm; see Wu et al. 2005), the hydrodynamical evolution should be Γ 2 ∝ R −m , and R ∝ t 1/(m+1) (Blandford & McKee 1976), with self-similar index m as where denotes the radiative efficiency of the shock, R the radius of the shock wave, Γ the bulk Lorentz factor of the shock, and t the observed time after the initial burst (Wu et al. 2005). The radiative efficiency can be considered as a constant at the early fast-cooling regime, and it evolves slowly after entering the slowcooling regime at later times. However, as long as we have power-law-distributed electrons, N e (γ e ) ∝ γ −p e , with γ e as the Lorentz factor of the electrons and N e the electron number, and the distribution index p does not deviate too much from 2, the effect of the evolving is not significant (Wu et al. 2005). Because the typical value of p for relativistic-shock-accelerated electrons is 2.2 (Achterberg et al. 2001, also see Curran et al. 2009 for observational constraints of an average value of k ∼ 2.24 with six Swift GRBs), which is close to 2, here we ignore the evolution of .
In this work, we focus on radio peaks occurring at T 90 + 0.5 days in the rest frame only, thus the effects of the equal-arrival-time surface (EATS) can be ignored. That is because, for a typical GRB, the radio-peak time t peak ∼ T 90 + 0.5 days < Γ −2 R/2c, which is the maximum time delay due to the EATS effects. Also, because the major concern of this paper is the behavior of GRB radio afterglows, we ignore the inverse Compton (IC) scattering, for IC is not a dominant factor at lower frequencies (Wu et al. 2005).
Following the procedures of Wu et al. (2005) and van der Horst (2007), we first consider the evolution of spherical shock waves with various k values. For synchrotron radiation caused by power-law-distributed electrons, a key time in the light curve is t cm , which marks the transition point between the early fast-cooling phase to the late slow-cooling phase, with an electron minimum Lorentz factor γ m equal to the electron cooling Lorentz factor γ c , and a minimum frequency ν m equal to the cooling frequency ν c (Sari et al. 1998;Panaitescu & Kumar 2000). We define ν cm = ν m (t cm ) = ν c (t cm ). The transition occurs at (2) where m e is the electron mass, m p the mass of a proton, e the charge of an electron, c the speed of light in vacuum, σ T the electron's Thomson scattering cross section, E cm the isotropic energy of the shock wave at t cm , and ζ 1/6 = 6 (p − 2) / (p − 1). Following the assumptions of Wu et al. (2005), we take typical values for parameters as follows: the isotropic energy E cm = 10 53 erg, ε e = 1/3, ε B = 10 −2.5 , and the electron distribution index p = 2.2. It should be noted that although some later works, including Nava et al. (2014), Beniamini & van der Horst (2017), as well as Santana et al. (2014), hinted that ε e ∼ 0.15 and ε B ∼ 10 −8 −10 −3 should be the case for some GRBs, the exact values of ε e and ε B do not change the overall shape of the afterglow light curves. Thus, our calculations still apply to their samples.
The characteristic flux density F ν,max at any given time can be obtained with the characteristic flux density of synchrotron radiation at a certain frequency on t cm as F ν,max (t cm ). (It should be noted that F ν,max is different from F peak,ν because the latter denotes the peak flux of the radio light curve at a given frequency ν, rather than the radiation peak across the radio spectrum). F ν,max can be calculated as follows, regardless of fast or slow cooling . (4) Also, at the radio band with low-enough frequencies, the synchrotron self-absorption (SSA) effects cannot be ignored. In this case, we have another characteristic time t a , which denotes the time when ν = ν a , where ν a is the self-absorption frequency. Because in this work the later-time behavior of the radio light curve is our main focus, we only discuss the later-time slow-cooling regime, with ν a = min {ν as,< , ν as,> }, and where c 0 ≈ 10.4 [(p + 2) / (p + 2/3)] ≈ 15.24, and the strength of the magnetic field B as well as shock radius R can be calculated with their corresponding values at t cm (e.g., see Equations (13) and (14) in (Wu et al. 2005)).
In the slow-cooling phase, the transition from ν as,< to ν as,> occurs at t am , when ν a = ν m . For various combinations of physical parameters, typically we have t am ∼ 10 3 t cm . For many cases, the GRB afterglow may have already finished the ultra-relativistic self-similar evolution and even began the transition to the Newtonian regime by the time of t am , making Blandford & McKee (1976) no longer applicable. And it is worth noting that in our numerical discussions of late-time afterglow, including the transrelativistic phase, in Section 2.1, the evolution of ν as,> is included in calculations, with corresponding effects of this factor in this regime considered. Thus, we do not take ν = ν as,> into consideration in analytical calculations.
For typical parameter values, we have ν cm ∼ 10 3 − 10 4 GHz, which is much higher than in the radio band; hence in this work, we only take the ν < ν cm band into consideration. Firstly, we calculate the light curves for radio afterglows. If the flux F ν ∝ t α has a positive temporal index α before a turning point and then shows a negative α, we define the turning point as a radio peak. We only consider the 0 k 2 case. The reason is that the self-similar solution proposed by Blandford & McKee (1976) only stands as long as k < 4. And in reality, it is difficult for the circumburst environment to show density distributions steeper than k > 2, since the original medium distribution and the stellar wind from the progenitor are the only main factors affecting k.
Thus, similar to Wu et al. (2005) and van der Horst (2007) by comparing the values of characteristic times t a , t cm , t c , as well as t m , with the latter two items corresponding to the time points when ν c and ν m equal to the observing frequency ν, respectively, we can get the overall profile of GRB radio afterglow light curves. Depending on the values of k, three types of light curves can be obtained. For smaller k with 0 < k < 32−12p−8 12−3p+(p−4) , when ν a (t cm ) < ν < ν cm , with ν a (t cm ) means the SSA frequency at t cm , we have t a < t cm < t m < t c , and the  Huang et al. (1999). Here we take the GRB isotropic energy Eiso = 10 53 erg, characteristic density A = 1 cm −3 , and adiabatic shock wave = 0. The y-axis are shown in arbitrary units. It can be seen that for the k = 0 case, the low-frequency afterglows peak during the transrelativistic-to-Newtonian phase, and the relations between peak times/fluxes and frequencies are different from the ultrarelativistic regime, while for the k = 1.1 and 2 cases, the radio afterglows never enter the trans-relativistic phase during our calculations. Figure 3. The upper panel shows the temporal evolution of Γv/c with various isotropic energies. The transrelativistic evolution featuring 0.2 < Γv/c < 2 is marked with the gray-shaded region bounded by the dotted lines. Here we take k = 0 and n = 1 cm −3 . It can be seen that the lower the Eiso is, the earlier the transrelativistic phase begins. The middle panel shows the temporal evolution of Γv/c for various circumburst densities, with k = 0, Eiso = 10 53 erg. It can be seen that the larger the density n, the earlier the transition begins. The lower panel shows the temporal evolution of Γv/c with various k, Eiso = 10 53 erg, and characteristic density A = 100 cm −3 . It can be seen that the smaller the k is, the earlier the transition time.
peak occurs at t m , with flux before and after peak as For ν < ν a (t cm ), the order changes to t cm < t a < t m < t c , and the flux peaks again at t m . The light curve around peak time can be expressed as It should be noted that, in this case, the upper limit of k, i.e., 32−12p−8 12−3p+(p−4) , depends on the electron distribution index p, as well as the shock radiation efficiency , and is quite sensitive to both parameters. For p ∼ 2.2, ∼ 0, this upper limit is ∼ 1.037. When the value of k increases to 32−12p−8 12−3p+(p−4) < k < 4/3, at higher frequencies, ν a (t cm ) < ν < ν cm , we have t a < t cm < t m , with t c disappearing in the radio band. In this case, the high-frequency light curve still peaks at t m , and the radio flux around this time can be calculated as (9) For ν < ν a (t cm ), t cm < t a < t m , t m remains as t peak,ν , and the peak-time light curve is (10) And both the high-and low-frequency radio fluxes peak at t m again.
For a larger k (4/3 < k < 2), the SSA frequency at t ac is lower than ν cm , that is, ν ac < ν cm , and the radio band should be divided into three sections. For ν ac < ν < ν cm , t a < t c < t cm < t m , while we have t c < t a < t cm < t m for ν < ν a (t cm ). In these two cases, the slope between t cm and t m is slightly negative as long as > 0, only with an early peak at t a . However, a significant break does exist at t m , as can be seen from the equation below. Considering possible observational errors, such a turnover could still be taken as peak time: (11) While for ν < ν a (t cm ), we expect t c < t cm < t a < t m , again with t m serving as peak time. Here the light curve can be expressed as (12) It can be seen that the GRB radio afterglow peaks at t peak = t m , regardless of k value or observing frequency. Because we have t m ∝ ν −2(m+1)/(4m+k) , the peak times and peak fluxes can be calculated with Eqs. 7-12 as That is, the a and b indices we seek should be For = 0, it can be seen that t peak,ν ∝ ν −2/3 . Even for > 0 cases, as long as the shock radiation efficiency is low enough, the t peak,ν ∝ ν −a index a should not deviate too much from 2/3. For the F peak,ν ∝ ν b index b with = 0, b increases with k, that is, b = 0 when k = 0, and b = 1/3 when k = 2. If > 0, the b value is larger than that in the = 0 case, although such a difference is not significant for larger k. In Figure 1, the relations between k and b with various are shown, according to which the k value can be drawn from b. And it should be noted that for later peaks occurring at T 90 +several days, = 0 usually applies, and thus a single k value can be constrained with b, as shown in Eq. 16.
Eqs 15 and 16 cannot be applied for high-frequency (ν ac < ν < ν cm ) radio afterglows with k > 4/3, the flux of which peaks at t a 1 day. Because for real GRBs radio follow-up observations usually do not begin at such early times, and the behaviors of afterglows are often contaminated by the early reverse shocks, such a high-frequency early peak with a large k have little value for our purpose, and its behavior is beyond the scope of this work.

Late-time Afterglow Behaviors and Transition from Relativistic to Newtonian Phase
The GRB afterglow begins the transition from relativistic to Newtonian phase (the so-called "transrelativistic" stage) when the shock wave expands to a larger radius. The evolutionary phase with Lorentz factor 0.2 Γv/c 2 can be considered as "transrelativistic", because by this time, the shock velocity has been significantly decreased and the amount of swept-up circumburst medium becomes large enough, while the afterglow evolution in this stage can be described by neither Figure 4. The temporal evolution of characteristic frequencies for various k (lower panel). Here the dashed, dotted and dashdotted lines represent the evolution of cooling frequency νc, minimum frequency νm, as well as SSA frequency νa. The solid lines show the evolution of multi-band radio peak frequencies, as well as peak times. It can be seen that considering trans-relativistic effects and evolution of νa, the peak time t peak,ν for these cases should be ≈ ta. However, other possibilities for peak also exist for various physical parameters.
the relativistic self-similar solution proposed by Blandford & McKee (1976), nor the Newtonian Sedov-Taylor solution (Taylor 1950;Sedov 1969). Here, numerical calculations based upon dynamical equations proposed by Huang et al. (1999), which provide a unified framework from relativistic to deep Newtonian stages of shock-wave evolution, are required, with examples shown in Figure  2. And because in late phases GRB afterglows should be adiabatic, it is safe to assume = 0 in our transrelativistic calculations. Figure 3 shows the temporal evolution of Γv/c for adiabatic shocks exhibiting various physical parameters, with the trans-relativistic phase occurring during T 90 +dozens to hundreds of days. And considering the transrelativistic phase only begins when the swept-up mass is comparable with the mass of the shock wave itself, it is clearly seen that a larger E iso , a smaller A, and a larger k can lead to a later transition.
According to numerical calculations, for spherical ejecta, the radio flux peak no longer appears at t m during such late times, due to the change of blast-wave dynamics, as well as the transition from ν a = ν as,< to ν a = ν as,> , which may occur before or after the be-ginning of transrelativistic phase, depending on the circumburst density profiles and E iso . For example, as can be seen from Figure 4, t a acts as the peak time for our typical parameters, though other possibilities also exist. Hence, the values of a and b in this phase differ significantly with predictions made by the self-similar shock-wave theory.
In late times covering the transrelativistic regime, the analytical expressions to describe the evolution of the bulk Lorentz factor Γ and other quantities no longer exist. However, the relations between t peak,ν , F peak,ν and ν can be analyzed numerically. We fitted a and b with various parameter values and found that both a and b show no strong dependence on E iso or A. While these indices do change more significantly with k. For typical parameter values, the larger the k value, the flatter the late time a, and the smaller the b. Generally speaking, for shocks with E iso = 10 53 erg, A = 100cm −3 , = 0, p = 2.2, ε e = 1/3, and ε B = 10 −2.5 , we have a ∼ 1.6, and b ∼ 0.95 for the k = 0 case. While a decreases from ∼ 1.6 to ∼ 1.0, b decreases from ∼ 1.0 to ∼ 0.89 when k increases from 0 to 2. By comparing Fig. 4   3(c), it can be seen that the late peaks due to the transition of ν a could appear earlier than the beginning of the transrelativistic phase for the circumburst environment with k 1 and small A. However, because our numerical calculations cover the complete from as early as < T 90 + 0.1 days, the general trend mentioned above still applies.
We also investigate the effects on transrelativistic a and b caused by the electron distribution index p, the proportion of electron and magnetic field energies ε e and ε B , as well as shock radiation efficiency . We find that for electrons with distribution index 2.0 < p < 2.5, the values of a and b depend weakly on p, and a larger p can lead up to a steeper b, as well as a shallower a. However, such a trend is not as clear as the k dependence mentioned above, and the change in a is hardly noticeable. For example, a decreases from 1.6 to 1.55 when p increases from 2.1 to 2.4. Meanwhile, a larger p can give rise to a b value as large as ∼ 1.1. Besides, both a and b depend weakly on ε B , and a smaller ε B can lead to flatter a and b. For example, a decreases from ∼ 1.6 to ∼ 1.4 when ε B decreases from 0.01 to 0.001, while neither ε e nor shows a strong influence on a and b. In short, the variation of b within a wide parameter range is relatively small, while the value of a is largely determined by k. With the value of a in the transrelativistic regime, one can put a stringent constraint on k.
For wideband observations covering ∼ 1−10 2 GHz frequency range, the values of the t peak,ν −ν and F peak,ν −ν indices a and b fitted from peak data should fall between a ∼ 2/3, b ∼ 0 − 1/3 for the ultrarelativistic regime, and a ∼ 1 − 0.9 as well as b ∼ 1.6 − 1 for the transrelativistic regimes, depending on k ∼ 0 − 2. Similar to Eqs. 15 and 16, such a transrelativistic trend with steeper a and b compared with the ultrarelativistic case can also be used to put constraints on k qualitatively, using a and b derived with multiband light-curve peaks.

Corrections for the Jet
Many GRBs have already entered postjet break evolution when reaching radio afterglow peak times, especially for low observing frequencies, that is, t peak,ν > t b , where t b denotes the jet break time (Chandra & Frail 2012). It is necessary to discuss the effects brought by a GRB jet in order to get a full understanding of multiband afterglows. Assuming the GRB jet-opening angle Figure 6. Theoretical GRB radio afterglow light curves considering forward-reverse shock interactions with Eiso = 10 52 erg(left) and 10 53 (right) erg, respectively. The other parameters are adopted as the typical values listed in Section 2. It can be seen that such interactions mainly affect higher frequency light curves at early times ( 1 day), with a smaller Eiso can give rise to an earlier FS/RS peak. Also, For frequencies < 4.8 GHz, the FS/RS peak is hardly noticeable, with the only prominent structure as the later peak shown in Fig. 2. is θ j , with an initial opening angle θ 0 . When the bulk Lorentz factor of GRB shock wave decreases to Γ < θ −1 j , the jet-opening angle becomes larger than the relativistic beaming angle 1/Γ, and thus, approximately speaking, the flux density from the jet is ∼ Γ 2 θ 2 0 times smaller than that from spherical shock.
Firstly, we consider a jet without sideways expansions, that is, θ j = θ 0 . Here we perform a similar analysis to Mészáros & Rees (1999). Assuming the flux density from spherical shock evolves as F ν,1 ∝ t a1 , a jet with the same parameters can give rise to a flux density F ν,2 ∝ t a2 . Because in the ultrarelativistic phase the bulk Lorentz factor evolves as Γ ∝ t (3−k)/(8−2k) , and F ν,1 /F ν,2 ∝ t a1−a2 ∝ Γ −2 θ −2 0 , it can be seen that That is, analytically speaking, for a relativistic jet without sideways expansions, the slope of the afterglow light curves during t > t j is the same as the slopes presented by Eqs. 7-12 minus a correction factor 3−k 4−k . It can be seen that no radio peak exists after t j , no matter the k value. If a jet exists during the radio afterglow evolution, observationally the radio-peak time t peak,ν should equal the jet break time t j , if no peak occurs before t j , and t peak,ν no longer depends on the observing frequency ν in this case, making a = 0. If we assume that t j < t m , and that t j is not far from t m , we can calculate that b = 1/3. Such a b value cannot be distinguished from peaks caused by spherical shocks. Although Huang et al. (2000) show that the shock dynamics can be affected by the jet, according to the numerical calculation, the afterglow flux decreases smoothly after jet break when sideways expansions of the jet can be ignored even with evolving temporal index, and later emission peaks do not appear. Thus, the analytical conclusions on a and b still apply.
Analytical expressions for jets with sideways expansions only exist when k = 0, according to Sari et al. (1999). In this case, we also have a 2 < a 1 , and a 2 0. That is, no emission peak exists after the jet break. Hence, similar to jets without expansions, the possible peak marks the jet break time t j , with a = 0, b = 1/3 (providing that t j occurs not too early). For arbitrary k values, a numerical calculation based upon dynamical equations provided by Huang et al. (2000) is needed. We take the local speed of sound c s = c/ √ 3 as the speed of jet expansion to calculate the dynamical evolution and find that the existence of a jet has a greater impact on late-time afterglow behaviors. During this time, the condition Γ > θ −1 j can be satisfied once more, making the afterglow flux rise again, forming a new emission peak. As early peaks are produced by spherical shocks, the late peaks due to jet effects appear earlier at higher frequencies, too. Because for typical parameters the late peak appears at several hundred days after GRB prompt emission, the analytical solutions proposed by Sari et al. (1999) or Rhoads (1999) do not apply, even for the k = 0 ISM environment. The flux level of such a late peak is much lower than that of early peaks, for the afterglow has already entered transrelativistic phase during this time, making the relativistic beaming effect much less significant. And it should be noted that when k ∼ 2, we can almost have Γ > θ −1 j satisfied all the way. As a result, it is possible that no jet break exists in the stellar-wind environment. Figure 5 shows the temporal evolution of the jet-opening angle for various k, and its relation with 1/Γ. According to numerical results, the a value of the late emission peaks for sideways-expanding jets decreases with increasing k, which is similar to the case with spherical shock waves. The value of a decreases from ∼ 1.1 to ∼ 0.88 when k increases from 0 to 2. Meanwhile, the value of b for jets increases with larger k, which is opposite to the spherical case. And b is quite small this time, even smaller than the ultrarelativistic case. For example, our results show that b ∼ 0.15 for k = 0, and b ∼ 0.45 for k = 2. However, it should be noted that all such results are based upon one single assumption, that is, the jet expands at a constant speed of c/ √ 3. In reality, this assumption may not be satisfied after the afterglow enters the transrelativistic regime. Latetime afterglows may expand quite slowly, making the jet-opening angle θ j < 1/Γ, thus no late peak can be formed. Also, earlier numerical simulation shows that jets do not have significant sideways expansions . In one word, effects caused by jet may be more complex in real GRBs.

Sample Selection
We collected data on GRBs with multiband radio light-curve observations (that is, GRBs with radio light curves in no fewer than two bands) from the literature, in order to put constraints on the circumburst density profiles based on Eqs. 13 and 14. Only the observed flux and corresponding epochs, rather than fitted peak values, such as the ones shown in Table 4 of Chandra & Frail (2012) or Table 1 of Li et al. (2015), were utilized in this analysis. Because it can be seen from Eqs. 13 and 14, as well as in Section 2.1 that theoretically speaking, both relativistic and transrelativistic cases give rise to nonnegative a and b, no matter of the k value, high-frequency GRB radio afterglow light curves should exhibit a higher F peak , as well as a t peak no later than lower-frequency peaks, a general trend of later and weaker peaks at lower frequencies is required during our sample selection, in order to comply with the theoretical framework of this analysis. For this reason, GRB 980519 (Frail et al. 2000b) was excluded, due to the late appearance of higher-frequency peaks compared with the lower-frequency ones, which is hard to explain with the GRB afterglow model based upon synchrotron radiation, unless we take large observation errors and light-curve modulation due to interstellar scintillation into consideration. Besides, GRBs 050416A (Soderberg et al. 2007), 020813, and 050713B  were omitted, due to their higher peak fluxes compared with their lower-frequency ones. Although our strategy is somewhat biased, such a bias could be considered as originating from model constraints, especially for our radio-peak-focused approach to get the k value.    Notea "VLA" for Very Large Array; "JVLA" for the upgraded Karl G. Jansky Very Large Array; "ATCA" for Australia Telescope Compact Array; "GMRT" for Giant Meterwave Radio Telescope; "WSRT" for Westerbork Synthesis Radio Telescope; "Ryle" for Ryle Telescope; "AMI" for Arcminute Mircrokelvin Imager. b Error bars of t peak,ν denotes the times on which the adjacent measurements are taken around the corresponding peak of the light curve. c Flux decreases monotonically during the complete observing campaign without a definite peak; we took the highest Sν as peak flux. d Flux increases monotonically during the complete observing campaign without a definite peak; we took the highest Sν as peak flux. e The 4.9 MHz light curve of GRB 980329 exhibits two peaks, according to Taylor et al. (1998). Since the first peak appears earlier than the 8.3 MHz peak, which is in contradiction with the standard afterglow model, here we adopt the later peak for our analysis. f The 30 GHz data acquired by the Owens Valley Radio Observatory's 40 m telescope from Galama et al. (2000) were discarded due to scintillated light curve without a definite peak, and the later appearance of the highest flux compared to the 15 GHz data; the 1.43 GHz data acquired by VLA from Galama et al. (2000) were discarded due to the earlier peak time compared to higher frequencies. g The 22.46 GHz data acquired by VLA from Berger et al. (2001) were discarded due monotonic light curve and later appearance of the highest flux compared to the 15 GHz data. h The 8.46 GHz data point with Sν = 566 ± 34 µJy acquired at 2000 Octpber 5.216 (UTC) by VLA from Harrison et al. (2001) was adopted as the peak flux; a later rebrightening at October 12.771 with Sν = 644 ± 126 µJy was discarded due to the later appearance compared to the 4.86 GHz peak, as well as the larger error bar. i The 1.43 GHz data acquired by VLA from Berger et al. (2003) were discarded due to scintillated light curves and earlier appearance of the highest flux measurement compared to the 4.86 GHz data; the 22.5 GHz data were discarded due to monotonic light curve and a maximum flux significantly lower than the 8.46 GHz data. j Two flux measurements with similar flux levels, Sν = 728 ± 55 and 749 ± 63 µJy, were performed at 2004 January 4.33 and January 15.35 (UTC), respectively, as shown in Soderberg et al. (2004b). It seems that the radio afterglow of this burst suffered significantly from interstellar scintillation, according to data presented in Soderberg et al. (2004b), thus making the exact peak time somewhat hard to determine. k The 15 GHz measurements acquired by the Ryle Telescope from Soderberg et al. (2006) were discarded due to large measurement errors and lack of data points; the 22.5 GHz measurements acquired by VLA from Soderberg et al. (2006) were discarded due to lack of data points. l An 8.46 GHz measurement performed at 2007 October 12.04 (UTC) shows a flux level of Sν = 431 ± 51 µJy (Perley et al. 2008), which is slightly larger than the peak data listed in this table. The earlier Sν = 430 ± 50 µJy was chosen as a peak considering the theoretical predictions for earlier peaks at higher frequencies. However, due to no 8.46 GHz data has been taken between the two measurements, it is difficult to know the exact peak time for the GRB 071003 radio afterglow at this frequency. m ACTA 5.5 GHz data for GRB 100418A listed in Moin et al. (2013) were discarded, due to the monotonic light-curve behavior, which is inconsistent with the VLA data.  2015) were discarded due to the lower flux and later peak compared to the 4.8 GHz data. p The 13 and 15 GHz measurements acquired by JVLA from Cucchiara et al. (2015) were discarded due to the lower flux compared to the 7 GHz data. q The 1.255 GHz measurements acquired by JVLA from Laskar et al. (2016) were discarded due to a peak flux that appeared earlier than the 1.644 GHz peak, with a flux higher than the 5.0 GHz peak. r The 1.45 GHz measurements acquired by JVLA from Alexander et al. (2017) were discarded due to an earlier peak compared to 1.77 GHz; all data with observed frequency > 13.5 GHz were discarded due to monotonic decline with maximum Sν lower than 11.0 GHz peak, due to the lack of early observations.
In total, we selected 32 GRBs, including 31 long bursts and 1 short one (GRB 170817A), as our sample. For observations at certain frequencies in our sample, GRBs do not comply with the theoretical trend of higher F peak , earlier t peak at higher ν -such data points have been omitted during our analysis. Meanwhile, we largely ig-nored the light-curve peaks occurring earlier than ∼ 1−2 days in the observer's frame to avoid possible contamination from early reverse shocks, because as can be seen in Fig. 6, the typical peak time due to forward-reverseshock interactions should be at ∼ T 90 + several×10 −1 days, which corresponds to ∼ T 90 + 1 − 2 days in the observer's frame because the typical redshift value of GRB is z ∼ 3 (Jakobsson et al. 2005). Also, 50 GHz data are not taken into consideration, as in such high frequencies, the t peak = t m conclusion does not apply (see Section 2). Considering the sparsely sampled nature of GRB radio observations, as well as the not-sosteep slopes of temporal evolution for radio afterglow light-curve fluxes (see Fig. 2), we took an approximate approach to make our analysis: The maximum observed flux is recognized as of the peak flux F peak,ν at a certain observing frequency, the corresponding observing time as the estimated peak time t peak,ν , and the epochs of adjacent observations give the upper and lower limits of t peak,ν . Besides, in order to make constraints with as many observing frequencies as possible, we took another strategy to utilize the monotonically evolved light curves: the reading of the last data point of a monotonically increasing light curve was taken as the lower limit of F peak , and the second-to-last epoch as the possible lower range of error for t peak . Similarly, for a constantly declining light curve, the first data point sets the lower limit of peak flux, and the second observing epoch the upper limit of peak time.
Some GRBs have been observed by various telescopes at similar frequencies. For example, GRB 030329 has been observed by WSRT at 1.4 GHz (van der Horst et al. (2005), and GMRT at 1.28 GHz (van der Horst et al. (2008)). Also, GRB 100418A has been observed by VLA and ACTA at 4.95 and 5.5 GHz, respectively Moin et al. (2013). Because calibrations between different instruments could lead to extra complications, in this case, we select peak data for analysis from as few telescopes as possible. Hence, our main uncertainties arise from errors in t peak,ν estimation, especially for GRB light curves with fewer observed data points. Information on the t peak,ν and F peak,ν for each band of our GRB samples, as well as descriptions of the discarded outlier data, is listed in Table 1. It is worth noting that as seen in Fig. 3, considering the long-lasting process of transition between relativistic to deep Newtonian shock waves, nearly all of our samples should be still in the relativistic or transrelativistic phase when the light-curve peaks occur.

Constraints on Circumburst Environment
The values of the a and b indices of 32 sample GRBs are listed in Table 2. Figure 8 shows the a values, observing frequencies, as well as peak times in light curves of each sample GRB. Figure 9 shows the corresponding b value, the observing bands, along with peak fluxes. And a statistic of the k value in the circumburst environment for our complete sample can be found in Figure  7. We compare these fitted indices with theoretical predictions from Eqs. 15 and 16 and Section 2.1, taking the circumburst density data provided by Chandra & Frail (2012) as well as possible observational errors into consideration. It can be seen that more than 40% of our samples can be explained with uniform ISM distribution, while nearly one-fifth samples can be described with stellar-wind environment. And due to various observational complications, both k = 0 and k = 2, or any 0 < k < 2 circumburst density distributions work for about 30% of the samples, while two GRBs cannot be properly constrained with our analysis using multiband radio peaks.

GRBs with ISM Density Profiles
According to our method, a total number of 14 GRBs out of the 32 samples can be solely explained with a k = 0 ISM environment. As listed in Table 2, we find that the values of a and b for multiband radio afterglow peaks of GRBs 980329, 980703, 000301C, 020903, 031203, 070125, 100418A, and 171010A lie between the ultra-and transrelativistic cases, and thus can be explained with transrelativistic afterglows under an ISM environment. And as can be seen in Figs 8 and 9, although due to observational uncertainties, the time-frequency index a shows large fitting errors for GRBs980329, 020903, and 100418A, which may lead to multiple possible k values, the flux-ν index b can be constrained stringently enough. In fact, as noted in Section 2.1, considering the late occurrence of the transrelativis-tic phase for k = 1.1 or 2 cases, a b ∼ 1 value alone can still lead to a confident conclusion of k = 0, no matter of the a value. Besides, tt is worth noting that the radio peaks of GRBs 020903 and 070125 show "breaks" of a and b values at a certain frequency, which could imply a transition from ultra-to transrelativistic phases, as the high-frequency peaks could appear during the ultrarelativistic phase with flatter a and b values, while the lowfrequency ones occur later with steeper transrelativistic indices.
Compared with existing works using multiband modeling, it should be noted that Berger et al. (2000), Soderberg et al. (2004b), and Chandra et al. (2008) have drawn similar conclusions for GRBs 000301C, 031203, and 070125, respectively, that is, all these three GRBs occur in ISM environments. And Starling et al. (2008) Figure 7. A statistic on the k distributions for the 32 GRB sample. It can be seen that the radio peaks of 14 GRBs can be explained by the k = 0 density profile (with six samples showing possible contamination by the reverse shock). Another six samples exhibit behaviors compatible with predictions for the k ∼ 2 case. Also, we have 10 GRBs that can be explained by density distributions between the ISM and stellar-wind environment, along with another 2 bursts occurring in environment hard to constrain, due to observation errors.
found a possible ISM explanation for GRB 980329 and demonstrated that GRB 980703 can be explained by both k = 0 and 2 environments. Frail et al. (2003) also showed that the wind model cannot be ruled out (although not favored) to fit the light curve of GRB 980703, which is still compatible with our constraints. An exception is GRB 171010A, which has been explained by Bright et al. (2019) using a steep density profile, although some hard-to-explain features do exist. This conclusion is somewhat inconsistent with our constraints, which could be due to the relatively large fitting errors for the b index, thus leaving large room for alternative explanations.
However, although the rest six GRBs in this category, including GRBs 050820A, 051022, 130427A, 141121A, 160509A, and 160625B generally follow predictions of the k = 0 ISM environment, Cenko et al. 2006, Perley et al. (2014, Cucchiara et al. (2015), along with Alexander et al. (2017) all pointed out a possible reverse-shock contribution of the radio afterglow peaks in these bursts. With the existence of reverse-shock contamination, the multiband peak method of analysis may no longer apply, thus the reliability of our related results could be compromized.

GRBs Compatible with a Stellar-Wind Environment
GRBs occurring in k = 2 density profiles as implied by multiband radio peaks include GRBs 980425, 000418, 011121, 030329, and 100814A, along with the only short GRB in our sample, GRB 170817A. Among these, the behavior of the GRBs 980425 and 030329 light-curve peaks follows predictions for late-time afterglow (ν a = ν a,> ) in the wind medium. And Berger et al. (2003) and van der Horst et al. (2005) have fitted the afterglow light curves of GRB 030329 with both ISM and wind models, whose results are compatible with our conclusions.
Meanwhile, the b indices of GRBs 000418 and 031203 are relatively shallow and are consistent with expanding jets in the ISM environment. However, their a indices are too small to be explained with jets. It also should be noted that for several bursts, the medium-and lowfrequency radio peaks may be explained with jet breaks, although the details depend much on the jet-opening angles (hence the jet break times). Analysis with multi-band radio peaks may only hint at multiple explanations, and data from other bands (e.g., optical) should be utilized to clarify the existence of jets. Besides, it is quite possible that the speed of the jet's sideways expansion during the transrelativistic phase should be slowed down, making the results drawn in Section 2.2 no longer stand.
The value of a index of GRB 011121 is consistent with the existence of a jet, while the b index falls within predictions for the k = 2 case. Considering the simultaneous appearance of the multiband peak and the large errors in peak-time estimations that can be attributed to observational effects, a burst occurring in a wind density profile is a more suitable conclusion, which is supported by Price et al. (2002c). A similar case applies for GRB 100814A, for which De Pasquale et al. (2012) has pinned down to the existence of a jet.
For the gravitational-wave-associated GRB 170817A, the only short burst in our sample, its multiband radio afterglow exhibits a shallow a with a shallow b. Such a case is more consistence with relativistic jet, which is compatible with the off-axis jet identified by previous works (e.g., see Fong et al. 2019).

GRBs Explained by k between 0 and 2
This category consists of 10 samples, which are GRBs 970508, 970828, 990510, 991208, 000926, 010921, 021004, 060218, 111215, and 140304A. The uncertainties here mainly arise from the relatively large errors for peak-time/flux estimates. For example, the peak behavior of GRBs 970508, 990510, and 010921 lies somewhere in between the k = 0 trans-relativistic and k = 2 late time (ν a = ν a,> ) cases, when taking fitting errors into consideration, making any 0 k 2 possible. This is compatible with the conclusion drawn by Starling et al. (2008), who found a homogeneous circumburst medium for GRBs 970508 and 990510 with X-ray, optical, and IR afterglows. However, the exact k for each sample is hard to pin down precisely, due to large uncertainties in a and b.
For GRBs with high circumburst densities, including GRBs 991208, 000926, 021004, 060218, and 111215A, the indices of a and b can be explained with later transrelativistic shocks in the k ∼ 0 environment. However, because a larger A can lead to a later t cm (thus a later peak time t m in radio afterglow), the larger k possibility, which can produce a steeper b during the relativistic phase, cannot be excluded for these bursts. When compared with existing works, Harrison et al. (2001) showed that the temporal behavior of GRB 000926's light curve can be better modeled with a Comptonized uniform medium; Soderberg et al. (2006) fitted GRB 060218 with both k = 2 and k = 0 density profiles, as well as low an E iso ; and Zauderer et al. (2013) has pointed out that the circumburst medium of GRB 111215A should be in the form of stellar wind, which are all consistent with our analysis.
For GRBs 970828 and 140304A, although the a and b indices can be explained with both k = 0 and 2 possibilities when taking fitting errors into consideration, it should be noted that the afterglows of these bursts were significantly affected by reverse shock, as shown by Djorgovski et al. (2001) and Laskar et al. (2018). And the radio light curve of GRB 140304A might even be contaminated by interstellar scintillation (Laskar et al. 2018). All of these factors can complicate our analysis, and thus compromise the robustness of the results.

GRBs Hard to Explain with Multiband Light-curve Peaks
Finally, we have two outliers that cannot be well constrained by the multiband peak method. One of them is GRB 000911, whose a and b values are much steeper than theoretical predictions for both trans-or ultra-relativistic shock waves, thus making its afterglow behavior hard to explain, even considering the fitting errors. Another hard-to-explain case is GRB 071003, which is located near the Galactic plane, with its radio data severely affected by scintillation (Perley et al. 2008), thus leaving a large room for errors in the a and b indices, therefore making multiple explanations, including both relativistic or transrelativistic ejecta propagating in the k = 0 and k = 2 circumburst environment, possible.

SUMMARY AND DISCUSSION
In this paper, we investigated the radio afterglows of gamma-ray bursts occurring in a power-law-distributed circumburst medium with n ∝ R −k , and an arbitrary k. We show that one can use multiband radio afterglow peak time t peak,ν and F peak,ν data to put a constraint on the density-distribution index k. We find that in the relativistic phase of the afterglow evolution, the peak time t peak,ν corresponds to t m , that is, the time when the minimum frequency of electrons equals to the observing frequency. And we have t peak,ν ∝ ν −a , a ∼ 2/3, which is independent of the shock radiation efficiency, while for F peak,ν ∝ ν b , the value of b depends more strongly on k. For adiabatic shocks, b increases from 0 to 1/3, if we change k from 0 to 2. And in the transrelativistic phase, similar dependencies between peak time/flux and frequency exist, with steeper a and b values than in the ultrarelativistic phase. Figure 8. The relations between radio afterglow peak fluxes t peak,ν and observing frequencies ν/νmin for our sample GRBs. The dashed line with data points in each figure shows the fitting results from observations, the blue dotted lines denote the possible range of fitting errors, while the black solid lines from top to bottom represent the theoretical predictions for relativistic jet break, relativistic spherical shock wave, k = 2 late time (νa = νa,>), and k = 0 transrelativistic cases, respectively. . The relations between radio afterglow peak fluxes F peak,ν and observing frequencies ν/νmin for our sample GRBs. The blue dashed line with data points in each figure shows the fitting results from observations, the blue dotted lines denote the possible range of fitting errors, while the black solid lines from top to bottom represent the theoretical predictions for transrelativistic, k = 2, = 0 relativistic, and k = 0, = 0 relativistic cases, respectively. It should be noted that as seen in Section 2.1, the b value for transrelativistic cases is insensitive to k.
It can be seen that by comparing multiband radio-peak data, one can determine the value of the circumburst medium distribution index k. We carry out an analysis for 32 GRBs with multiband radio afterglow light-curve data, and find that half of them can be explained with transrelativistic afterglows under uniform ISM environments. Only ∼ 20% of our sample GRBs can be certainly determined that they occur in wind-like environments with larger k, although it should also be noted that the radio peaks of nearly 30% of our sample GRBs can be interpreted by a k value between 0 and 2. Most of these results are compatible with existing analyses of individual bursts. However, Yi et al. (2013) found that nearly all their 19 sample GRBs show signs of k ∼ 0.4 − 1.4, with a typical value of k = 1, by analyzing early forward-reverse-shock evolution. Obviously, our conclusion is different from that of Yi et al. (2013). One of the possible explanations could be that the emission regions of GRB afterglows become larger after entering the self-similar phase. And circumburst medium distributions at these regions may be quite different from the early forward-reverseshock-dominated areas, which are much closer to the GRB central engine. That is because the mass-loss processes of GRB progenitors induced by the stellar wind or other mechanisms have limited influence and may not change the density distribution farther away. Besides, the sample utilized by Yi et al. (2013) does not overlap with our radio sample. It is difficult to judge whether such inconsistency is due to diversity among individual bursts or changes in density distribution at different distances around a single GRB. If complete followup observations for one burst become available in the future, from early X-ray and optical detections to late-time radio observations, we can use these information to fully understand the GRB circumburst medium distributions at various distances and possible changes within.
As seen in Section 3, the main uncertainty in our analysis comes from errors in peak-time/flux estimation. As for each GRB we consider the maximum observed flux as its approximate peak flux and the corresponding observing time as an estimated peak time in our index fitting, while the real peak may be situated at any time between the two observing times adjacent to our estimation, with a flux no lower than our estimation, if the sampling of the light curve is not frequent enough during the follow-up observations, especially if the data points near the peak time are sparsely distributed, the fitted a and b values using our method do have larger errors, and the circumburst medium distribution cannot be pinned down without doubts. However, it is worth noting that although such uncertainty does exist, our constraints on k are still better than those from single-band light-curve slope fitting. On one hand, in many cases, it is nearly impossible to fit the light-curve slopes reliably in order to calculate k with only a handful of data points; while our approach described in Section 3 can show some preliminary results at least. On the other hand, even if we can make enough observations for one GRB in a single band, those data far away from the emission peak often exhibit large errors due to their low fluxes and instrument limitations, and a precise fitting on k cannot be guaranteed. On the contrary, with more data on hand, the uncertainties on peak estimation are greatly reduced this time. Hence, with more high-quality light-curve observations made available in the future, the parameters on multiband radio peaks can be better determined, and more stringent constraints on circumburst density distributions can be obtained.
Because multiband afterglow observations covering the low-frequency range may show signs of transrelativistic transitions, the values of a and b at various bands could be quite different. Because existing radio observations usually sample a few frequencies only, we can only see signs of such transitions in a handful of bursts. Thus, the uncertainties here cannot be ignored. If one can carry out radio follow-ups at more frequencies in the future, such transitions should be unveiled. Besides, the new generation of large radio telescopes, including the Five-hundred-meter Aperture Spherical radio Telescope (FAST) and the upcoming Square Kilometer Array (SKA), have much higher sensitivities at lower frequencies, and thus can be used to detect GRB late-time afterglows (e.g., see Li et al. 2015;Zhang et al. 2015;Ruggeri & Capozziello 2016), and sample more frequency bands at longer wavelengths in order to get a more complete picture of transrelativistic shock behaviors.