Prospects for Detecting γ-Ray Bursts at Very High Energies with the HADAR Experiment

Recent ground-based observations of TeV photons have significantly deepened our understanding of the nature of gamma-ray bursts (GRBs). However, many fundamental problems remain unsolved concerning the physical mechanisms behind GRBs, necessitating the need for sufficient statistical data. The High Altitude Detection of Astronomical Radiation (HADAR) experiment utilizes a wide-angle water Cherenkov telescope, presenting a novel approach to measure the spectra and variability of GRBs from 10 GeV to 10 TeV energy ranges with unprecedented photon statistics and thereby break new ground in elucidating the physics of GRBs, which is still poorly understood. In this study, a time-dependent numerical modeling technique is utilized to simulate extensive light curves and spectral energy distributions of synthetic GRB afterglow emissions. By harnessing the remarkable capabilities of HADAR, we evaluate its potential in detecting GRB afterglow emissions at energies >10 GeV. Through our analysis, we unveil the prospect of detecting an estimated 5.8 GRBs annually, facilitating a systematic investigation into their reliance on model parameters. Future HADAR observations would offer valuable insights into the magnetic field and the environmental conditions surrounding GRBs.


Introduction
A gamma-ray burst (GRB) jet running into an external medium is expected to generate luminous GeV-TeV emission lasting from minutes to several hours (see Zhang 2018 and references therein).The high-energy emission might result from different kinds of radiation mechanisms, such as the electron-synchrotron emission (Zhang & Mészáros 2001), the synchrotron self-Compton (SSC) emission of electrons (Wei & Lu 1998;Sari & Esin 2001), the proton-synchrotron emission (Böttcher & Dermer 1998;Totani 1998), and some other hadron-related emission (Pe'er & Waxman 2005).The highestenergy observation of GRBs holds a key for understanding the physical process of acceleration, energy transfer in collisionless astrophysical shocks, and so on.Besides GRB physics , such measurements will also provide crucial diagnostics of ultrahigh-energy cosmic-ray and neutrino production in GRBs, advance observational cosmology by probing the high-redshift extragalactic background light (Inoue et al. 2010) and intergalactic magnetic fields (Dzhatdoev et al. 2020), and contribute to fundamental physics by testing Lorentz invariance violation (Shao & Ma 2010;Acciari et al. 2020) with high precision.
Until now, seven GRBs have been reported with sub-TeV to TeV emissions, 11 which not only shed new light on GRB studies but also pose challenges.MAGIC observed a hardening of the spectrum of GRB 190114C, along with the dip measured in LAT flux.This provides evidence for the presence of an additional SSC radiation component besides synchrotron emission (MAGIC Collaboration et al. 2019a, 2019b).On the other hand, GRB 190829A significantly exceeds the synchrotron burn-off limit, challenging the simple synchrotron or synchrotron plus SSC model (H.E.S.S. Collaboration et al. 2021).This requires a more complex and elaborate model for explanation.The mechanism responsible for high-energy photon emission in GRB afterglows is still controversial, as they may consist of multiple components depending on the fireball parameters (Wei & Lu 1998).Therefore, having sufficient observations is crucial.
Furthermore, sub-TeV to TeV emissions in GRBs exhibit remarkable longevity during the afterglow phase, accompanied by unusually high energy photons.Notably, TeV photons from GRB 180720B and GRB 190829A persist for extended periods, spanning tens of hours and nearly 3 days after the initial burst (Abdalla et al. 2019;H.E.S.S. Collaboration et al. 2021).In a recent breakthrough, the LHAASO detector detected the highest-energy photons ever recorded, reaching an astonishing 18 TeV, during the observation of GRB 221009A (Huang et al. 2022).Just as impressive, the Carpet-2 air-shower array has also reported an air shower consistent with a 251 TeV photon energy 4536 s after the GBM trigger (Dzhappuev et al. 2022), although it remains uncertain whether the observed 251 TeV photon is from GRB 221009A or either of these Galactic sources (Sahu et al. 2023).To comprehensively understand the underlying physical processes and determine the prevalence of such high-energy emissions across all GRBs, it is imperative to accumulate a sufficiently large sample of events (Berti & Carosi 2022;Miceli & Nava 2022).
The complexity of recent observations on very-high-energy (VHE; >100 GeV) GRB scenarios highlights the need for ongoing and persistent observations to improve our understanding of GRB physics.The Fermi-LAT has made significant contributions to the study of GRBs in the MeV-GeV energy range since its operation (Ajello et al. 2019).Due to its effective area of about 7000 cm 2 at 1 GeV and decreasing effectiveness at higher energies, the detection of sub-TeV and higher-energy emissions is challenging for the LAT.In contrast, the effective collecting area of ground-based facilities increases with energy.Imaging atmospheric Cherenkov telescopes (IACTs) like MAGIC, H.E.S.S., and VERITAS have provided valuable observations in the highest energy range but are limited by their narrow field of view (FOV).The recent LHAASO experiment has demonstrated the advantages and capabilities of extensive air-shower arrays (EASAs) for observing GRBs.For a comprehensive comparison of key space and ground facilities in GRB science, please refer to Tsvetkova et al. (2022) for a review.It is crucial to have experiments with unprecedented photon statistics and sensitivity to investigate the parameter space of a wide range of VHE-transient emitters and their characteristics.More facilities, especially those with large FOV and high sensitivity, may provide important simultaneous observations for spectral energy distribution (SED) modeling.
The High Altitude Detection of Astronomical Radiation (HADAR) is this kind of experiment, harnessing the advantages of both IACTs and EASAs.In this study, we investigate HADAR's capability to observe GRBs and the information it provides about them.This paper is organized as follows: Section 2 provides a brief description of the HADAR array and an overview of the experiment's current status.In Section 3, we present the GRB afterglow model and describe our method for sampling the GRB population.Our results and conclusions are presented in Sections 4 and 5, respectively.

HADAR Experiment
The HADAR project seeks to establish an imaging air Cherenkov detector array with a wide FOV, specifically designed to detect gamma rays above 10 GeV.As an innovative observation technique, it incorporates four largedimensional hemispherical glass shells, filled with purified water, as light collectors to capture Cherenkov photons produced by gamma-ray showers in the atmosphere (Cai et al. 2017;Chen et al. 2019).
In order to validate the effectiveness of this new technique, a prototype water lens, with a diameter of 0.9 m, was installed and operated at the Yangbajing Observatory in 2015 and successfully detected cosmic-ray signals (Cai et al. 2017;Chen et al. 2019).Furthermore, the construction of a second prototype of a 2.5 m water lens commenced in 2021 and is scheduled to begin operation shortly.The primary scientific goal of HADAR is to cover γ-ray astronomy in the energy range from 10 GeV to 10 TeV, providing exceptional sensitivity while maintaining a broad FOV.Further details regarding HADARʼs instrument layout, configuration, and performance analysis under different incident rays can be found elsewhere (Xin et al. 2021).Our Monte Carlo simulations have yielded some notable outcomes, including the expectation of observing two to three prompt emission GRBs per year (Xin et al. 2021).

Model and Methodology
In this section, we will provide an introduction to the model used for GRB afterglow emission, as well as the distributions of model parameters that we will employ.We will then outline the methodology utilized to generate a diverse population of GRBs, with each GRB possessing its own distinctive light curve and SED.Subsequently, we will simulate the observation of these synthetic GRBs using the HADAR experiment, considering various observational conditions.This comprehensive analysis allows us to evaluate the sensitivity of HADAR across different scenarios and acquire valuable insights into its observational capabilities.

Time-dependent Numerical Modeling
We investigate an impulsive outflow characterized by its kinetic energy, E k,iso , and initial Lorentz factor, Γ 0 .This outflow propagates into a constant-density external medium, denoted as n 0 .The nonthermal electrons initially follow a power-law distribution with an index of p.These electrons lose energy through synchrotron emission in the presence of a turbulent magnetic field and through inverse Compton scattering of the self-produced synchrotron photons.The equipartition factors for the energy in electrons and magnetic field in the shock are denoted as ò e and ò B , respectively.The lower-energy emission primarily arises from the Thomson regime, with the Klein-Nishina effect considered in the inverse Compton process (Wang et al. 2010).
Since the hydrodynamic properties of the forward shock responsible for the afterglow have been extensively discussed in previous literature, including the formulation governing the initial electron distribution, photon production, and absorption processes (Sari et al. 1998;Panaitescu & Kumar 2000;Sari & Esin 2001), we refrain from restating them here.Interested readers can refer to those works for further details.In the Appendix, we provide a concise overview.It should be noted that we do not take into consideration the angular extent or any angular structure of the jet.Additionally, we assume no lateral spreading of the blast wave.While the reverse shock emission is expected to contribute to early GeV emission, we have not considered it in this study.
Following the standard afterglow model (Sari et al. 1998), the light curve at a given observed frequency (ν) can be calculated as F(t, ν) = F(t, ν, E k,iso , Γ 0 , p, ò e , ò B , n 0 ).To model the GRB afterglow emission and obtain the time-evolving SEDs, we utilize a numerical code developed by Liu et al. (2013).

GRB Population Sampling
To ensure dependable forecasts concerning the detection rates of GRBs with the HADAR experiment, we have produced a population of GRBs, each with its unique time evolution calculation.This entails employing a set of specific parameters (E k , Γ 0 , ò e , ò B , p, z, n 0 ).It is noteworthy that these parameters are not universally applicable and showcase a significant variation among the GRB samples.Thanks to the comprehensive efforts in establishing GRB afterglow parameters through earlier research, the remaining parameter space is manageable in this work.To produce the GRB population, we conducted an exhaustive literature search to provide values for the parameters.The key model parameters requested for afterglow emission of the GRB population are concisely explained below and are visually depicted in Figure A1.
(1) E k , the initial kinetic energy of the blast wave, is often linked to E iso through the relationship η = E iso /(E k + E iso ).E k = 5 E iso is assumed to determine the values of E iso (Beniamini et al. 2016;Li et al. 2023).The distributions of E iso are derived from measurements obtained by the Swift satellite, as documented in the GRB table.12(2) Γ 0 , the initial bulk Lorentz factor, has been estimated using various methods, including the opacity method (Abdo et al. 2009;Ackermann et al. 2014), the afterglow onset method (Molinari et al. 2007;Ghirlanda et al. 2012Ghirlanda et al. , 2018)), and the photosphere method (Mészáros & Rees 2000).Among these approaches, the afterglow onset method provides a more reliable estimation with fewer uncertainties.In our investigation, we adopt the distribution of Γ 0 derived from a robust sample of 66 GRBs with measured t peak as the standard distribution (Ghirlanda et al. 2018).Additionally, we utilize their upper limits on Γ 0 as a test distribution to assess the impact of this parameter.
(3) ò e and ò B , two microphysical parameters, have been extensively studied in various works.Santana et al. (2014) have presented compelling evidence revealing that ò e exhibits a remarkably narrow distribution.This distribution spans only one order of magnitude, ranging from 0.02 to 0.6, indicating that a minority of GRBs possess values of ò e below 0.1.In contrast, the distribution of ò B is considerably wider than that of ò e , lacking a distinct range.Santana et al. (2014) found that previous afterglow modeling studies primarily focused on the range of ò B between 10 −4.5 and 10 0 , whereas their own results derived from X-ray and optical observations indicate that the distribution of ò B is concentrated within the range from 10 −8 to 10 −3 .We adopt their X-ray-derived results as the standard distribution for ò B , while also considering the results from optical observations and multiple works by other researchers as test1 and test2 distributions.(4) Typically, the electron index p predominantly ranges from 2 to 3 (Shen et al. 2006;Wang et al. 2015).In this study, we utilize the inferred distribution of the electron spectral index p from Wang et al. (2015) with subsamples within the gold sample, as well as their best Gaussian fits, as the standard p distribution for our sampling.Furthermore, we select their full samples as the test p distribution for our analysis.
(5) Previous studies on afterglows have indicated that a constant-density medium is generally more appropriate (Yost et al. 2003;Schulze et al. 2011).In our analysis, we employ the same value of n 0 for each group of samples, with 1 cm −3 being the standard choice.However, to assess the impact on the observation rate, we vary this value and examine the cases of 0.1 and 30 cm −3 .(6) Wanderman & Piran (2010) studied the intrinsic redshift distribution; we employ their redshift distribution as the standard distribution.
High-energy photons from cosmological emitters can still suffer the γγ interactions with the extragalactic background light (EBL) before reaching the observer.The collective emission of any high-energy emitting cosmological population will exhibit an absorption feature at the highest energies.We consider several EBL models, including those proposed by Gilmore et al. (2012)

Simulating HADAR's Observations of GRBs
A database of GRBs, including their afterglow emission F(t, ν) across a broad range of energy and time spectra, is acquired through randomly sampling from the GRB population mentioned above.We are now poised to identify which one could be observed by HADAR.The positions of the GRB population are distributed uniformly in the sky, and the explosion times are also distributed uniformly over the course of 1 yr.The observation zenith angle for each burst by HADAR should be less than 30°, while simultaneously ensuring that the Sun and the Moon are below the horizon.
Then, we calculate HADAR's sensitivity flux for each GRB.This result is much more refined, accurate, and closer to reality compared to the previous approach of taking the average of all GRBs.For ground-based experiments, a 5σ deviation from the background distribution is required to claim the discovery of a GRB.A power-law distribution of background cosmic rays is assumed.For each pseudo-GRB, the signal significance is calculated by = S N N s b , where thus, the sensitivity is

. 3 sen
The performance of the detector, including its effective area (A eff ) at different energy ranges and observation zenith angles (θ), angular resolution (Ω), and proton-γ discrimination capability (η), has been adopted from a previous study (Xin et al. 2021).
The energy range considered in this study spans from 10 GeV to 10 TeV, and the maximum observation time after the burst is 10 days.A GRB is included in the observational list if its signal significance, through either time integration or differential time binning, exceeds 5σ.

Results
After conducting comprehensive simulations of the afterglow emission of GRBs, involving the generation light curves and SED of a large sample of thousands of GRBs, and simulating the observations of each GRB using the HADAR experiment, we have obtained valuable insights.Now, we can perform a statistical analysis on all the GRBs observed within the sensitivity range of HADAR.This analysis will allow us to quantitatively evaluate the observational capabilities of HADAR in detecting and studying GRBs.

GRB Detection Rate with HADAR
Figure 1 presents a visual representation of the simulated results obtained from the observations conducted during the HADAR experiment.The focus is specifically on the light curve and energy spectrum of a single GRB.In the left panel of Figure 1, it can be observed that the assumed location of the GRB falls within the FOV of HADAR, and at the assumed burst time, neither the Moon nor the Sun is present.However, approximately 1000 s after the burst, the conditions for HADAR observation become unfavorable.This could be due to either the burst's zenith angle exceeding 30°or the presence of the Moon or Sun.It is only on the second night that the burst position reenters the HADAR FOV.However, by this time, the decaying flux of the GRB has become too low to be effectively observed.The light curve of the GRB and HADAR's sensitivity to its observation over a span of 10 days are depicted in this figure.It is important to note that only the first day of data can be collected by HADAR.The right panel of Figure 1 displays the average spectrum of the GRB within the initial 1000 s, with the dashed line indicating the impact of EBL absorption.This figure serves as an illustrative example, and the same procedures are applied to all other synthetic GRBs by comparing their fluxes with HADAR's sensitivity.
According to the sensitivity limit of current space-born GRB detectors, an ideal imaginary 4π all-sky detector on average is roughly two to three GRBs per day (Zhang 2018).To estimate this rate conservatively, we employed Poisson sampling with a mean value of 600 GRBs per year.In total, we conducted 1000 simulation observations with HADAR.In the left panel of Figure 2 we present the distributions observed in our experiments, which were derived using the EBL model proposed by Gilmore et al. (2012).Additionally, we have included a Poisson fit with a mean value of 5.8.Among these, approximately 2.1 and 5.77 GRBs are observed within the first 1000 s and the first day, respectively.The results with different EBL models are presented in the right panel of Figure 2. Notably, under the EBL model proposed by Helgason & Kashlinsky (2012), the maximum annual GRB observations  reach 7.9.Conversely, the default EBL model used in our study yields the lowest prospects for GRB observations.

Systematic Uncertainties
To ensure the reliability of the calculations presented in Figure 2, we conducted further analysis on the GRB observation rate of HADAR across various parameter distributions.Specifically, we substituted one standard parameter distribution with a test distribution, as illustrated in Figure A1. Figure 3 presents the outcomes obtained when employing the test distribution.Notably, it can be observed that while variations in parameters such as Γ 0 and p have minimal impact on the annual observation rate, replacing the environmental medium density n 0 and magnetic field factor ò B does yield noticeable effects.When the test2 magnetic field factor is adopted, as shown in Figure A1, and a stronger magnetic field is present, it results in a reduction in the number of HADAR GRB observations.This can be easily understood, as a larger amount of energy is radiated into the synchrotron radiation compared to the SSC components.When the density of the medium is set at 1 cm −3 , the annual observation rate of GRBs by HADAR is approximately 5.8.Reducing the density by a factor of 10 leads to a decrease in the observation rate to ∼3.3 per year, while increasing the density by 30 results in an observation rate of 4.8 per year.When the external medium increases to 30 cm −3 , the reverse Compton scattering internal absorption becomes severe, resulting in a decrease in the flux of GRBs.As a result, the number of observed GRBs is actually lower compared to when the external medium is 1 cm −3 .These results underscore the potential of future HADAR observations in offering valuable insights into the shock magnetic field and environmental conditions surrounding GRBs.

Summary
Multiwavelength observations of GRBs, particularly with a focus on TeV energy range radiation, play a critical role in advancing our understanding of GRB characteristics.The HADAR instrument, employing a refractive water lens, offers unique advantages with its wide FOV and substantial effective area in the sub-tens of GeV energy range.In this study, we conducted numerical calculations to determine the detailed synchrotron and SSC emissions using a parameter set obtained from previous statistical studies.With an assumed annual allsky burst event of 600 GRBs, we estimated HADARʼs detection capacity, yielding an expected detection rate of approximately 5.8 GRBs per year with n 0 = 1 cm −3 .
To account for uncertainties arising from parameter distributions, we also explored alternative distributions for one of the parameters and compared the resulting flux with HADARʼs sensitivity using an extensive statistical sampling under specific conditions.Our findings highlight the potential of future HADAR observations in providing valuable insights into the shock magnetic field and the environmental factors influencing GRBs.Considering an impulsive relativistic outflow with an initial kinetic energy E k and initial Lorentz factor Γ 0 , propagating into an external medium with constant density n 0 , the overall dynamic evolution of radius (R), swept-up mass m, and the Lorentz factor Γ (Huang et al. 2000) can be described as respectively, where m p is the proton mass, c is the light speed, b = -G - 1 s s 2 is the velocity of the shock, and its relationship with the bulk Lorentz factor is given by Figure A1.Key parameters utilized in generating GRB populations, with red lines representing standard distributions and blue and dark-green lines representing test distributions used to evaluate differences resulting from varying parameter distributions.Top left: E iso , derived by using the redshift-observed data of 400 GRBs.The fluence measurements are obtained from the Swift-BAT instrument, and the data are sourced from the GRB table available at https://swift.gsfc.nasa.gov/archive/grb_table/fullview/.The observations included in the analysis extend up to GRB 230506C.Top right: distribution of Lorentz factor, as reported by Ghirlanda et al. (2018).Middle: distributions of microphysical parameters ò e and ò B , adapted from Santana et al. (2014).Bottom left: distribution of redshift, adopted from Wanderman & Piran (2010).Bottom right: distribution of electron index p (Wang et al. 2015).
with the adiabatic index ˆ( ) g = G + G 4 1 3 .Relativistic shocks serve as sites for particle acceleration, magnetic field amplification, and photon radiation (Zhang 2018).In this paper, we consider ò e and ò B as constants, representing the fractions of shock internal energy distributed to electrons and magnetic fields, respectively.The distribution of accelerated electrons is commonly modeled as a single power-law function of the electron Lorentz factor γ e , expressed as g g µ - dN d e e p , within the range of minimum and maximum Lorentz factors (Blandford & McKee 1976;Sari et al. 1998).The minimum Lorentz factor, γ m , can be described by the following equation: where p represents the electron index, σ T is the Thomson scattering cross section, and the magnetic field intensity (Sari & Esin 2001;Panaitescu 2005) is parameterized as


The radiative processes affect the shape of the electron distribution, resulting in a broken power-law spectrum with several segments (Sari et al. 1998;Sari & Esin 2001).The cooling Lorentz factor of electrons, γ c , due to synchrotron and SSC radiation (Sari et al. 1998;Sari & Esin 2001), is described as Additionally, the electron Lorentz factor is limited up to Here the parameter Y evaluates the effect of SSC cooling on the synchrotron spectrum (Sari & Esin 2001).The electrons lose energy through synchrotron emission in the presence of a turbulent magnetic field and through inverse Compton scattering of the self-produced synchrotron photons.Synchrotron radiation power is computed using the traditional (Crusius & Schlickeiser 1986) formula Here ( ) = + + -+ -+ F q g q q q q q , 2 ln 1 2 1 1 qg gq 4 2 1 4 2 .The intrinsic spectral flux in the observer frame can be derived from the radiation power using the following equation: In this equation, D L represents the luminosity distance between the source and the observer.The value of D L can be calculated using cosmological parameters, specifically Λ M = 0.27, Ω Λ = 0.73, and H 0 = 71 km s −1 Mpc −1 .
To summarize, the calculation of the light curve for an individual GRB involves considering a specific set of model parameters, namely E k , Γ 0 , ò e , ò B , p, z, n 0 , in relation to the observed frequency F(ν, t).
In Figure A1, we illustrate the distributions of the essential model parameters utilized for sampling the GRB population.For detailed descriptions of these parameter distributions, please refer to Section 3.2.

,
Inoue et al. (2013), Stecker et al. (2012), Finke et al. (2010), and Helgason & Kashlinsky (2012); the EBL attenuation introduced in the work of Gilmore et al. (2012) is used by default.The values of the optical depth in a designated energy range at a specific redshift are interpolated.

Figure 2 .
Figure 2. Right: the detected count of GRBs utilizing HADARʼs observation sensitivity, considering an intrinsic average annual burst of 600.GRBs are generated with standard parameter distributions.A total of 1000 simulated experiments were conducted, and a Gaussian fit is employed for analysis.Gilmore et al. (2012) is considered.Left: four different EBL models are considered; they are from the works of Inoue et al. (2013), Stecker et al. (2012), Finke et al. (2010), and Helgason & Kashlinsky (2012).

Figure 3 .
Figure3.The observed count of GRBs using HADAR is investigated, employing a methodology similar to that demonstrated in Figure2, to explore the impact of varying a specific parameter on the distribution.Each panel is accompanied by the Poisson fit results of the respective distribution.The red line corresponds to the result obtained in Figure2.In the top left panel, the standard distribution of Γ 0 is replaced with an alternative test distribution.Similarly, in the top right panel, the standard ò B distribution is varied by utilizing a test distribution.Moving to the bottom left panel, the standard p distribution is substituted with a test distribution.Lastly, the density mediums are represented by different colored regions in the bottom right panel, where red, blue, and green correspond to density values of 1, 0.1, and 10 cm −3 , respectively.
the usual synchrotron function containing the modified Bessel function.In the same way, the total SED of the SSC radiation is