Contribution of γ-Ray Burst Afterglow Emissions to the Isotropic Diffuse γ-Ray Background

The isotropic diffuse γ-ray background (IGRB) serves as a fundamental probe of the evolution of the extreme Universe. Although various astrophysical sources have been proposed as potential contributors to the IGRB, the dominant population is still under debate. γ-ray bursts (GRBs) are among candidate contributors of IGRB, although they are not as frequently discussed as blazars or starburst galaxies. Recent observations of TeV emission from GRB afterglows have provided fresh insights into this subject. This work aims to investigate the potential contribution of GRB afterglows to the IGRB under the standard afterglow model. We carefully examine the influence of each microphysical parameter of the afterglow model on this contribution, with a particular emphasis on the significant role played by the initial kinematic energy. To determine the energy and quantify the contribution of GRB afterglow to IGRBs, we utilize the observed GRB afterglow energy emissions from the Swift X-ray Telescope and Fermi Large Area Telescope instruments. Our calculations, considering the synchrotron self-Compton emission, suggest that GRB afterglows make up less than 10% of the IGRBs. To enhance the precision of our findings, it is crucial to further constrain these parameters by conducting additional ground-based observations of GRB afterglows.


Introduction
The extragalactic γ-ray background was first detected by the SAS-2 satellite (Fichtel et al. 1977) and subsequently confirmed by the EGRET experiment utilizing a more sophisticated analysis involving the subtraction of individual sources and the Galactic emission (Sreekumar et al. 1998).The Fermi Large Area Telescope (LAT) experiment expanded the energy range from 100 MeV to 820 GeV with an improved event selection and a more precise removal of the cosmic-ray (CR) background and diffuse Galactic foreground, yielding a high-quality spectrum of the isotropic diffuse γ-ray background (IGRB; Ackermann et al. 2015).The IGRB is a nearly isotropic and diffuse emission that emerges from the combination of real diffuse γ-rays resulting from hadronic interactions of CRs with cosmic background radiation (such as cosmic microwave background or extragalactic background light, EBL; Berezinsky & Smirnov 1975) and unresolved extragalactic point sources such as blazars (Stecker & Salamon 1996;Mauro et al. 2014), millisecond pulsars (Faucher-Giguère & Loeb 2010;Siegal-Gaskins et al. 2011), misaligned active galactic nuclei (Inoue 2011;Mauro et al. 2013), and star-forming galaxies (Ackermann et al. 2012;Chakraborty & Fields 2012).Exotic physics, such as dark matter annihilation or decay, has also been considered as possible contributors to the IGRB (Taylor & Silk 2003;Liu et al. 2017;Hütten et al. 2018).Despite these advances, the exact composition of the IGRB remains one of the main intriguing questions in γ-ray astrophysics.Further observations and analyses are necessary to better understand this enigmatic phenomenon.
γ-ray bursts (GRBs) are generally believed to be the most energetic explosions since the Big Bang, and their high γ-ray luminosity suggests that they may contribute to the IGRB (Totani 1999).The calculation of Hartmann et al. (2003) on the prompt emission below MeV energy showed that GRBs can only account for a fraction of the observed IGRB at the soft γray band.Later dedicated calculations of the MeV-TeV flux during both the prompt phase and afterglow also suggested GRBs to be subdominant contributors to IGRB (Casanova et al. 2007;Ando et al. 2008).On the other hand, using a phenomenological model, recent studies have estimated the very high-energy (VHE) emissions during the prompt phase, demonstrating that GRBs could potentially make a considerable contribution to the high-energy end of the IGRB spectrum measured by Fermi-LAT (Yao et al. 2020).
While the VHE γ-ray emission is just a potential possibility in the prompt emission of a GRB, it is measured in afterglows of several GRBs.The dawn of VHE observations of GRBs begin with MAGIC and H.E.S.S. detection of VHE emission from GRB 190114C and GRB 180720B, respectively (Abdalla et al. 2019;MAGIC Collaboration et al. 2019).Later, other GRBs exhibiting VHE radiation were observed (Blanch et al. 2020a(Blanch et al. , 2020b;;Acciari et al. 2021).The latest LHAASO collaboration observations of GRB 221009A even detected photons with energy above 10 TeV (Huang et al. 2022).This is coupled with the fact that the VHE component of GRBs can be extended for several tens of hours after the initial burst, exemplified by GRB 180720B and GRB 190829A (H.E.S.S. Collaboration et al. 2021).Theoretical modeling has shown that the synchrotron self-Compton (SSC) mechanism (e.g., Fraija et al. 2019;Wang et al. 2019;Zhang et al. 2020;Joshi & Razzaque 2021), which results from the scattering of lowenergy photons by high-energy electrons, is a promising explanation for the VHE emission.Operation of the SSC mechanism does not require extreme conditions of a GRB and the circumburst environment, implying that the VHE component could be commonly present in GRB afterglows, although the absorption due to EBL and the limited sensitivity of instruments at the VHE band hinder a common detection among reported GRBs.Nevertheless, those undetected VHE photons of GRB afterglows including the reemitted emission by electron/positron pairs generated upon the absorption of VHE photons in the EBL are counted in the IGRB.These observational and theoretical results motivate us to further explore the contributions of GRB afterglows to the IGRB.
This paper is structured as follows.Section 2 presents the external shock model and radiation mechanism employed in this study.The sampling scheme utilized is described in Section 3.1.The calculation results are presented in Section 4, and the conclusions drawn from our analysis are summarized in Section 5.

Model Description
We employ a forward shock model utilizing a fireball at its center.This section encompasses detailed descriptions of the deceleration and radiation procedures involved.

Hydrodynamic Evolution of the Blast Wave of the GRB
We consider an outflow with an initial kinetic energy E k and initial bulk Lorentz factor Γ 0 that propagates into an external medium with a proposed constant density n for simplicity.The hydrodynamic process can be described by a differential equation, as presented by Huang et al. (2000): is the rest mass of the external medium for a constant density medium swept by the blast wave, and ò is a parameter describing the energy loss.The initial rest mass of the blast wave is denoted by The dynamic evolution of the outflow can be divided into three phases throughout the deceleration phase: the coasting phase, the Blandford-McKee self-similar phase, and the Sedov-Taylor self-similar phase.Equation (1) can describe all of these phases.To obtain the evolution of the Lorentz factor Γ of the blast wave, we should add the kinematics equation, as presented by Huang et al. (2000): is the velocity of the shock, and its relationship with the bulk Lorentz factor of the ejecta Γ is given by Blandford & McKee (1976): The adiabatic index, briefly denoted by G , equals 5/3 in the nonrelativistic case and 4/3 in the ultrarelativistic case.For each GRB afterglow, we model the evolution up to 1 yr after the explosion, which is sufficiently long to account for the VHE emission, given that the measured VHE afterglow lightcurve declines more rapidly than t −1 .

Shock-accelerated Electrons
In astrophysical systems, free particles undergo continuous acceleration, and their energy distribution can be described by the Fokker-Planck equation (Chandrasekhar 1943).In the context of GRB problems, this equation can be simplified as where N(γ, t) is the number density of electrons with energy γ at time t, ( )  g g is the rate of energy loss of electrons with energy γ, and Q(γ, t) is the injection rate of electrons with energy γ.
For a power-law distribution of electrons, the distribution of accelerated electrons in the slow cooling case is given by Sari et al. (1998): , .
In the fast cooling case, the distribution of accelerated electrons is given by , ., where R is the radius of the shock and n is the density of the ambient medium.

Afterglow Radiation Mechanisms
During the afterglow phase of GRBs, the outflow is subject to dynamic equations and accelerates electrons, which subsequently decelerate and radiate in the magnetic field of the shock.Assuming a strong shock and that a fraction of the energy transforms into magnetic energy with efficiency ò B , the jump condition yields the magnetic field energy density in the comoving frame e B B 8 , while e¢ is the energy density of the shock in the comoving frame, and it can be given by Blandford & McKee (1976), where n is the density of the ambient medium, m p is the proton mass, c is the speed of light, Γ is the bulk Lorentz factor of the outflow, and ĝ is the adiabatic index of the gas.
In the synchrotron radiation process, the minimum, cooling, and maximum characteristic Lorentz factors of electrons are given by where m e is the electron mass, t is the observer time, and σ T is the Thomson scattering cross section.The fraction of energy that the blast wave converts into electrons of the external medium is denoted by ò e .Additionally, Y ˜is the Compton parameter that reflects the ratio of the intensity of inverse-Compton radiation and synchrotron radiation, taking into account Klein-Nishina suppression.Understanding the minimum, cooling, and maximum characteristic Lorentz factors of electrons is crucial for interpreting the observed data and inferring the physical conditions of the GRB environment.After obtaining the electron spectrum at each time, we may calculate the synchrotron radiation and SSC radiation following the treatment detailed in Blumenthal (1970) and obtain the lightcurve of the afterglow of each simulated GRB.For simplicity, we ignore the equal-arrival-time effect.Adopting a flat ΛCDM Universe with H 0 = 71 km s −1 , Ω m = 0.27, and Ω Λ = 0.73, we can calculate the time-averaged spectrum of the simulated GRBs with a total number of 15,000 over a period of 30 yr.
To summarize the model description, there are seven parameters involved, as outlined in Table 1.With such a set of parameters, we may model the evolution of the GRB afterglow's radiation as a function of time, employing the code developed by Liu et al. (2013).In the subsequent sections, we will elaborate on the method of generating the GRB sample.

GRB Sample
To quantitatively estimate the overall flux contribution of GRBs to the IGRB, we generate a GRB sample following the observational properties of GRBs.Each GRB in this sample is defined by a specific set of physical parameters.Their lightcurves and spectral energy distributions are calculated numerically using the models described in Section 2. In this section, we begin by providing an overview of the method generating the GRB sample.We will also compare the statistical distribution of the radiation properties of the generated GRBs with that of the observed GRBs, so that we may verify the model parameters taken to generate the sample.

Sampling Scheme
The flux generated by the external-forward shock of a GRB jet is related to several parameters, including E k , ò e , ò B , Γ 0 , p, n, and redshift z.In order to obtain and reproduce a realistic GRB sample, taking into account their inherent variability, it is essential to possess knowledge of the distribution of each parameter.Fortunately, there have been previous statistical studies conducted on most of these parameters.In the subsequent sections, we will introduce the parameter distribution we have utilized for population sampling.

The Initial Kinetic Energy
To determine the initial kinetic energy of a phenomenon, we rely on its correlation with the prompt emission energy.This relationship is typically connected to the radiation efficiency of GRBs, which lies in the range of η γ = E γ,iso / (E γ,iso + E k,iso ) ∼ 0.01-0.9(Zhang et al. 2007).In our study, we assume that E k,iso = kE γ,iso , which means a radiation efficiency To obtain the initial value of E k,iso , we randomly sample a value from the distribution of E γ,iso .This distribution is represented by a Gaussian function that fits the logarithmic distribution of the GRB prompt emission energy from the Fermi Gamma-Ray Burst Monitor (GBM) catalog (Poolakkil et al. 2021), as depicted in Figure 1.A comprehensive overview regarding all Gaussian fittings related to the parameter distribution can be found in Table 1.The value of k (i.e., the radiation efficiency η γ ) is slightly varied to ensure consistency with the experimental energy distribution, and further details can be found in Appendix A. Normally, we accept k = 9.5 corresponding to a radiation efficiency η ∼ 10%.A maximum value of E γ,iso is set at 5 × 10 54 erg to avoid overly energetic GRBs, since the most energetic one is GRB 160625B with 5 × 10 54 erg by Band fitting (Poolakkil et al. 2021).We use a conservative approach to avoid the case that the most energetic GRB with the lowest redshift was coincidentally produced, which means that this one can serve as hundreds or thousands of them.

Energy Partition Factors
The parameters ò B and ò e are crucial for assessing the energy equipartition within the shock, representing the fractions of internal energy allocated to magnetic fields and electrons, respectively.Santana et al. (2014) has summarized and demonstrated that ò e exhibits a narrow distribution, varying by only 1 order of magnitude from 0.02 to 0.6, with few GRBs reported to have ò e < 0.1.In contrast, the distribution of ò B is significantly wider compared to ò e and lacks a clear range.
On the other hand, recent observations of GRBs with VHE emission have yielded smaller values of ò B compared to 10 −2 (Wang et al. 2019;LHAASO Collaboration et al. 2023).These findings align with theories that suggest SSC occurs when ò B < ò e (Sari & Esin 2001).In order to focus our analysis, we have adopted an ò B distribution derived from X-ray observations (Santana et al. 2014).Utilizing a Gaussian function, we have fitted this range to further sample the parameters.These parameter samplings are represented by the black lines depicted in the middle right panel of Figure 1, illustrating the distribution of ò B .To better understand the influence of ò B on the contribution to the IGRB, we have also tested the ò B distribution summarized from Santana et al. (2014), which follows a bimodal distribution and is represented by the red histogram in the middle right panel of Figure 1.For ò e , we use a Gaussian function to fit the distribution reported in Santana et al. (2014), which is also shown in the middle left panel of Figure 1.Regarding the test distribution of ò e , it is selected as the outcome of fitting the Gaussian distribution with the broadest possible broadening, which is also shown as the red line.When sampling ò e and ò B , we require that their sum, ò e + ò B , is less than 1.For the sake of simplicity, we make the assumption that these two microphysical parameters remain consistent and do not evolve with time.

The Initial Bulk Lorentz Factor and Medium Density
During the coasting phase, the ejecta maintains a constant bulk Lorentz factor Γ 0 before deceleration due to interaction with the external medium.Γ 0 represents the maximum value achieved by the outflow during this dynamic evolution and is a crucial yet challenging parameter to directly measure in GRBs.Various methods, such as the opacity method (Abdo et al. 2009;Ackermann et al. 2014), afterglow onset method (Ghirlanda et al. 2012(Ghirlanda et al. , 2018)), and photosphere method (Mészáros & Rees 2000), have been employed to estimate Γ 0 .Among these methods, the afterglow onset method offers a more reliable estimation with fewer uncertainties.In this study, we adopt the distribution of Lorentz factors Γ 0 estimated from a robust sample of 66 GRBs using the afterglow onset method (Ghirlanda et al. 2018) as the standard distribution.Additionally, we present the distribution of upper limits for 106 GRBs from the same study as a test distribution, represented by the red line on the top right of Figure 1.Throughout this paper, we assume an interstellar medium of a density of n = 0.1 cm −3 for the circumburst environment, on which the afterglow onset method relies.We also show the result for a higher density of interstellar medium n = 1 cm −3 in Figure 5 for comparison.

Index of Electron Injection and Redshift
The electron spectral index (p) for relativistic shocks resulting from first-order Fermi acceleration typically falls within the range of 2.2-2.4 (Achterberg et al. 2001;Spitkovsky 2008).Wang et al. (2015) provided the distribution of p for 85 GRBs and demonstrated consistent results among different subsamples.Among their gold sample of 45 GRBs, the distribution of p follows a Gaussian distribution with a mean of 2.33 and, on the other hand, with a mean of 2.58 for GRBs the X-ray and the optical emission of which is in the same spectral regime.We adopt the former distribution as the standard one for sample generation and take the latter one to test the influence of p, as shown in Figure 1.It should be kept in mind that we consider the sampling results with p > 2, which is normal in most GRB afterglow cases.For determining the cosmological distribution of GRBs, we rely on the intrinsic redshift distribution derived from the study conducted by Wanderman & Piran (2010).To make a comparison, we also show the observed redshift distribution calculated by the empirical broken power law from Yao et al. (2020), who had thoroughly cross-verified the redshift with observations, represented by the dashed line in Figure 1.

Consistency with Observed Energy of Afterglow
In order to validate our generated GRB sample, we compare the distribution of the isotropic-equivalent energy with that of the observed GRB sample.The latter consists of 31 Swift X-ray Telescope (XRT) GRBs, from which we calculate the isotropic γ-ray energy in the XRT band (0.3-10 keV) as a representation of the synchrotron band, and 11 Fermi-LAT GRBs with precalculated isotropic energy in the LAT band (0.1-10 GeV) hardly representing the SSC band, denoted as E XRT and E LAT , respectively.
During the afterglow phase of a GRB, the γ-ray flux increases sharply until the jet begins to decelerate.The deceleration time is denoted by t p and is typically at t dec E 150 10 erg 54 1 3 After the jet deceleration, the γ-ray lightcurve decreases with time more rapidly than t −1 , as observed in the afterglows of many GRBs, such as GRB 190114C, GRB 221009A, etc.Therefore, most of the γ-ray energy radiated in a GRB afterglow is released around the deceleration time.However, the onset of GRB afterglow is not always caught by instruments.We choose an onset time of 300 s for those GRBs without early observational data and calculate the total isotropic energy radiated in the afterglow by extrapolating the observed lightcurves to this time.As shown in Figure 2, we calculate the distribution of the ratio of energy released from 300 s, denoted by E γ (300 s), to that released in the entire afterglow phase E γ,tot for our synthetic populations of GRBs with the standard parameter distributions.We see that the values of E γ (300 s)/E γ,tot in about 80% of simulated GRBs are higher than 0.5, implying that E γ (300 s) could generally represent E γ,tot .We can take GRB 150314A as an example, which has a temporal index of −1.08 for XRT lightcurve.The isotropic-equivalent energy calculated is 2.2 × 10 52 (2.4 × 10 52 ) erg, corresponding to a 300 (200) s onset time.There is only a deviation of about 0.04 in log-scale axis, while the width of each bin is 0.75 in the left panel of Figure 3.This validates our choice of 300 s as the starting point to which we extrapolate the observed lightcurves and estimate the total energy released in the afterglow phase.For the case where t p > 300 s, the onset time is chosen based on the t p time.For the end of integral time, we accept the observation end time if the power of the lightcurve exceeds −1 to avoid divergent integration.The deviation from our estimation is minor, as predicted by the external shock model.
Figure 3 demonstrates that the energy distribution of the synthetic-generated GRBs, obtained by sampling with the aforementioned standard distributions of parameters, is well aligned with that of the observed samples.

Results
After obtaining a set of parameters for each GRB in the sample, we are able to calculate their multiwavelength afterglow lightcurves f (E, t).We then can get the time-integrated γ-ray flux of the afterglows in the generated samples.This will enable us to assess their overall contribution to the IGRBs.

Contribution from Afterglow
The total number of generated GRB afterglow data sets, denoted as N GRB , is estimated to be 1.5 × 10 4 .This corresponds to the number of GRBs observed over an approximately 30 yr period, assuming an average observation rate of 1-2 day -1 Figure 2. The energy release process during burst for our synthetic populations of GRBs with the standard parameter distribution.The y-axis represents the normalized cumulative distribution.The green and red solid lines denote the LAT band (0.1-10 GeV) and the XRT band (0.3-10 keV), respectively.We can see that most of the GRB population has a similar value close to 1 for E γ (300 s)/E γ,tot , which means E γ (300 s) is a good indicator of E γ,tot .
across the entire sky.To calculate the time-averaged spectral energy distribution, we employ the following equation: Here, ΔT is set to 10 9 s, providing sufficient duration to account for the time-decay flux of GRBs.The term e E in ( ) t g describes the internal absorption by ambient photons from the source.And high-energy photons will be attenuated and interact with EBL as they traverse intergalactic space.So the observed flux is more depressed than the averaged spectrum by . We consider two EBL models, including those proposed by Gilmore et al. (2012) and Inoue et al. (2013).These models offer valuable insights into the attenuation of high-energy photons during their journey to Earth.The absorbed high-energy photons may initiate electromagnetic cascades, which deposit the energy of these high-energy γ-ray photons to the GeV band.For the calculation of the electromagnetic cascade, we utilize the code developed by Liu et al. (2016).
Utilizing GRB samples generated with standard distributions, as depicted in Figure 1 and listed in Table 1, we show the time-averaged all-sky diffusive γ-ray intensity expected from GRB afterglows in Figure 4 considering η ∼ 10%, in comparison with the measured IGRB spectrum from Fermi-LAT.Absorption of the γ-ray flux by the EBL is important above several tens of GeV.The two EBL models employed do not yield different results, so we use the model by Gilmore et al. (2012) as the default EBL model for subsequent analysis.We also tried to employ a higher density of n = 1 cm −3 for the circumburst medium in Figure 5. Notably, upon considering EBL absorption, it is evident that the average flux contributed by the GRB samples reaches 6.5% (6.7% considering cascade) of the Fermi-LAT observed intensity of the IGRB at 10 GeV as an insufficient contributor.If we accept the SSC emission as an extra component and with the same definition as in Yao et al. (2020), R ext is 14% and β ext is −1.7 in this work without considering cascade.This implies that the afterglow contribution can be compared to the prompt phase in several GeV bands, it still holds a nondominant position.Considering suggestions that blazars (Ackermann et al. 2016) or starforming galaxies (Roth et al. 2021) dominate the IGRB, GRB afterglows may occupy a considerable fraction of the rest part of the IGRB.If the major component of the IGRB from blazars or star-forming galaxies can be accurately determined, the rest part of the IGRB may then put useful constraints on the properties of the GRB afterglow in turn.

Effect of Different Parameter Distributions
To assess the robustness of the results presented in Section 4.1, we examine the influence of various parameter distributions on the contribution of GRB afterglows to the IGRB.In the following analysis, we replace one of the standard   parameter distributions with the corresponding test distribution described in Section 3.1 one by one.Throughout the sampling procedure, we ensure that the GRB sample adheres to the observed distributions of E XRT and E LAT .This investigation allows us to evaluate the reliability of the results obtained.
In the first panel of Figure 6, the results obtained by employing the test distribution of Γ 0 (as shown in Figure 1) are presented.It can be observed that as the Lorentz factor transitions from the calculation using the exact transformed value of t p to higher values based on the upper limits calculated with optical emission, the time-averaged flux shifts slightly toward higher energy bands.For a higher initial Lorentz factor of outflow, it spends a little time to decrease sharply to a lower value, which hardly influences the time-averaged spectrum in the long term.
In the production of high-energy photons in GRBs, ò B plays a crucial role.The results of replacing the distribution of ò B are presented in the top right panel of Figure 6.The histogram sampling from the test distribution yields a much lower value of k = 1.2 (η ∼ 50%) to match the observed distribution of E XRT .Notably, the energy range of the IGRB we consider closely aligns with the typical energy of SSC emission.Noting the Compton parameter Y  (Sari & Esin 2001), it is evident that a higher value of ò B inhibits SSC emission to some extent.An increasing magnetic field equipartition factor leads to a softer spectrum and a weaker intensity of high-energy afterglow emission.But the lack of GeV afterglow weakly constrains the distribution of E k,iso and ò B , which introduces the most uncertainty for the results.The bottom left panel of Figure 6 presents the result obtained by replacing ò e , showcasing minimal variations in contributions derived from both the standard and test parameter distributions.The test distribution of index electron p implies k = 13 to match the energy constraints in Section 3.2, so it enhances the intensity from afterglow to some extent.
In this work, we do not consider the fluctuations induced by another distribution of the initial kinetic energy of the blast wave, which has been constrained by observations through the k factor in Section 3.2.However, we need to emphasize that the contribution from afterglow has a close relationship with the initial kinetic energy replaced by the k factor and high-energy population of it.If we add the short GRB in the GBM catalog (Poolakkil et al. 2021) for energy fitting, the distribution will become wider with σ = 0.92, which also increases the fraction of the higher energetic population and enhances the contribution.But for short GRBs, they always have a lower energy than the mean value of the Gaussian fit.So we choose to ignore the short GRBs to avoid potential errors.The high-energy cutoff in the distribution of E k,iso introduces fluctuations of about 1% into the results shown in Figure 4.It is evident that the entire distribution of the initial kinetic energy in logarithmic scale significantly impacts our results.If we have additional afterglow observations with long duration in the LAT band, we can obtain more compelling constraints for several parameter distributions, such as ò e and ò B .

Conclusion
It has been established that GRB afterglows are GeV-TeV γray emitters.In this work, we investigated the contribution of GRB afterglows to the IGRB with the standard forward shock model of GRB afterglows.To estimate this contribution, we synthesized a population of GRBs based on the energy budget distribution and spatial distribution derived from observed samples.A set of model parameters relevant for the afterglow emission is attributed to each GRB following the statistical analyses of these parameters proposed in previous literature.Through a detailed modeling of the lightcurve of each generated GRB afterglow, we are able to get the time-averaged spectrum and subsequently calculate the contribution to the IGRB.The dependence of the result on model parameters is also studied.Assuming η ∼ 10%, our result suggested that GRB afterglow contributes about 6.5% of the measured IGRB at ∼10 GeV.If the majority of the IGRB above 10 GeV arises from blazars or star-forming galaxies, as suggested by previous studies, the rest part of the IGRB actually provides a restrictive upper limit of the contribution from GRB afterglows, especially considering additional contributions, possibly from hadronic interaction of ultra-high-energy and exotic scenarios such as annihilating dark matter.If the rest part of the IGRB can be accurately separated from the majority, we may use it to constrain the general properties of some microscopic physical quantities of GRB afterglow among the population, which would benefit the study of ultrarelativistic shocks.

Acknowledgments
This work is supported by the National Natural Science Foundation of China (Nos.12075258, 12275279, U2031105) and the China Postdoctoral Science Foundation (No. 2023M730423)

Appendix A How to Determine the k Factor
We collect the isotropic energy emission for 31 Swift-XRT GRBs and 11 Fermi-LAT GRBs in their energy bands and draw two distributions in Figure 3.It is noted that we only choose the E LAT for those GRBs with long-term observation.In fact, we will have a precise sample to conform the E LAT distribution and a realistic contribution to the IGRB in the LAT band.Unfortunately, the long-term observations of LAT extending to thousands seconds toward GRBs are scarce.
Considering that E XRT has larger statistics and a more realistic distribution, we always demand the synthetic GRB samples according to the E XRT distribution.The width of distribution of the E XRT of the synthetic GRB sample is relevant to the width of E k after calculating, which demands E k to be as narrow as possible for a Gaussian: E log 52.94 0.7.
The peak of distribution of E XRT corresponds to the value of k.For ò B , a higher value will relatively facilitate synchrotron cooling and weaken the SSC process.For the standard parameter distribution we stated earlier, the value of k is 9.5 to fit the observed energy distribution, which is hardly changed if we alter the distribution of one of the parameters, except for ò B and the index of electron p.To show the contrast made by fluctuation of the ò B distribution, we accept another distribution, which is shown as red line in the middle right panel of Figure 1.This distribution corresponds to a lower value of k of 1.2 compared to the standard ò B distribution but still barely fits the XRT and LAT energy observations simultaneously.For a higher value density of medium n = 1 cm −3 , the radiation is increased for more electrons being involved in this process, corresponding to a lower value of k = 4.5.The test distribution of p demands k = 13.0, for a softer electron spectrum needs a lower radiative efficiency to fit the observation.In other calculations, k always remains 9.5.

Appendix B XRT and LAT Observed GRB Afterglow Energy Calculation
The total isotropic-equivalent emission energy in XRT band and LAT band in GRB afterglow in Figure 3. See Table B1 for more details.Notes.In this table, we show the X-ray energy emission in XRT band 0.3-10 keV of 31 GRBs and the γ-ray energy emission in LAT band 0.1-10 GeV of 11 GRBs in the afterglow phase.We use a single power law (SPL) to fit the lightcurve and calculate the isotropic-equivalent emission energy.For nine LAT GRBs, we choose them from the second LAT GRB catalog with a long observation term (Ajello et al. 2019).
m and γ c are the minimum and cooling characteristic Lorentz factors of electrons, respectively, and max g is the maximum Lorentz factor of electrons.The normalization factor is determined to satisfy the conservation of particle number,

Figure 1 .
Figure1.Key parameters we used for GRB sampling; the black (red) lines are the standard (test) distribution.Top left: histogram of energy distribution with Gaussian error is shown as black line; data are fromPoolakkil et al. (2021).Top right: the distribution of the Lorentz factor reported fromGhirlanda et al. (2018).Middle left: the ò e distribution.Middle right: the ò B distribution.Bottom left: the distribution of the electron index p.Bottom right: the distribution of redshift denoted as a solid line, adopted fromWanderman & Piran (2010).

Figure 3 .
Figure3.The distribution of γ-ray energy emission in the afterglow phase of 31 Swift-XRT GRBs (0.3-10 keV) and 11 Fermi-LAT GRBs (0.1-10 GeV) as well as our GRB sample with standard parameter sampling.Details on the calculation of experiment data can be seen in TableB1of Appendix B.

Figure 4 .
Figure 4. Comparison of calculated γ-ray spectra with the observed IGRB spectra from Fermi-LAT (Ackermann et al. 2015).The calculated flux is derived from GRB samples generated with standard parameter distributions.Long-dashed and solid lines represent spectra without and with EBL absorption, respectively.The EBL models of Gilmore et al. (2012) and Inoue et al. (2013) are considered.The increasing flux for secondary particles in the electromagnetic cascade, as predicted by the model proposed by Liu et al. (2016), is represented by a red solid line.The time-averaged spectrum ignoring the attenuation of the initial kinetic energy distribution without a high-energy cut is shown by a short-dashed line.

Figure 5 .
Figure 5.The time-averaged spectral energy distribution for different densities of circumstance.The wider solid lines represent the situation considering Gilmore's EBL model.

Figure 6 .
Figure 6.Comparison of calculated γ-ray spectra between the standard parameter distribution with one parameter replaced by the test distribution.Experimental data points are adopted from Ackermann et al. (2015).Dashed and solid lines represent spectra without and with EBL absorption, respectively.Top left: replacement of the standard distribution with the test distribution of Γ 0 .Top right: replacement of the standard distribution with the test distribution of ò B .Middle left: standard distribution with the test distribution of ò e .Middle right: standard distribution with the test distribution of p.

Table 1
Standard and Test Parameter Distributions a E k,iso = kE γ,iso .b Gaus (mean, sigma) are used to fit a histogram with a Gaussian function.c p > 2.