Jet Structure and Burst Environment of GRB 221009A

We conducted a comprehensive investigation of the brightest-of-all-time GRB 221009A, using new insights from very high-energy (VHE) observations from LHAASO and a complete multiwavelength afterglow data set. Through data fitting, we imposed constraints on the jet structure, radiation mechanisms, and burst environment of GRB 221009A. Our findings reveal a structured jet morphology characterized by a core+wing configuration. A smooth transition of energy within the jet takes place between the core and wing, but with a discontinuity in the bulk Lorentz factor. The jet structure differs from both the case of the short GRB 170817A and the results of numerical simulations for long-duration bursts. The VHE emission can be explained by the forward shock synchrotron self-Compton radiation of the core component, but requiring a distinctive transition of the burst environment from uniform to wind-like, suggesting the presence of complex pre-burst mass ejection processes. The low-energy multiwavelength afterglow is mainly governed by the synchrotron radiation from the forward and reverse shocks of the wing component. Our analysis indicates a magnetization factor of 5 for the wing component. Additionally, by comparing the forward shock parameters of the core and wing components, we find a potential correlation between the electron acceleration efficiency and both the Lorentz factor of the shock and the magnetic field equipartition factor. We discuss the significance of our findings, potential interpretations, and remaining issues.


INTRODUCTION
Gamma-ray bursts (GRBs) are among the brightest stellar-level events in the universe.The multi-wavelength afterglow radiation of GRBs spans a wide range of the electromagnetic spectrum, making them valuable subjects and tools in the study of time-domain astronomy.However, the understanding of the radiation characteristics and origins of GRBs and their afterglows in the very high energy (VHE; > 100 GeV) regime has been limited for a significant period of time because of lack of observations (Zhang & Mészáros 2001;Zou et al. 2009;Piron 2016;Nava 2018).After the launch of the Fermi detector, the Large Area Telescope (LAT) recorded a group of GRBs with photon energies ∼ 10 GeV up to 100 GeV, which opened a detection window in that energy range (Nava 2018).Additionally, detection of the VHE counterpart of GRBs has been a crucial scientific objective and challenge for imaging atmospheric Cherenkov telescopes (IACTs) for a long time (Acciari et al. 2011;Aliu et al. 2014;Alfaro et al. 2017;Berti & Group 2017;Hoischen et al. 2017).As a result, both observation and theory have been driving efforts to expand observations into the VHE window.
The Large High Altitude Air Shower Observatory (LHAASO) recently reported their observations of the VHE counterpart of GRB 221009A, claiming the detection of photons with energies exceeding 10 TeV (Huang et al. 2022).This represents the first instance of humans observing the VHE counterpart of a GRB in the energy range exceeding 10 TeV (Cao et al. 2023).Additionally, LHAASO is also the first detector to have monitored the VHE afterglow light curve of a GRB with both extremely high time and energy resolutions.LHAASO reported their observations and analytical conclusions in the 0.2 − 7 TeV energy range (LHAASO Collaboration et al. 2023).Thanks to the extremely high brightness of GRB 221009A, a wealth of afterglow data was collected through multiwavelength observations (e.g., Bright et al. 2023;Kann et al. 2023;Laskar et al. 2023;LHAASO Collaboration et al. 2023;Ren et al. 2023).Combining a rich assortment of multiwavelength afterglow data together with valuable insights obtained from the newly accessible VHE energy regime, the study of GRB 221009A has the potential to significantly enhance our understanding of GRB phenomena.Moreover, it will establish important connections with research domains associated with, for example, the physics of plasma under extreme conditions, the stellar formation and evolution, the ultimate phases of massive star lifecycles, and the environments in which they exist (Miceli & Nava 2022).
In this paper, using a dataset with both long duration time scale and wide energy range, we conducted an analysis of the radiation mechanism, jet structure, and burst environment of GRB 221009A.This paper is organized as follows.
In Section 2, we introduce our model considerations.In Section 4, we introduce our fitting parameters and method in detail.The results of model parameter inference and corresponding discussion are presented in Section 5. We summarize and conclusion the significance of our results in Section 6.We take the cosmology parameters as H 0 = 67.8km s −1 Mpc −1 , and Ω M = 0.308 (Planck Collaboration et al. 2016).

MODEL ASSUMPTION
Recent LHAASO observations of GRB 221009A have shown the presence of a remarkably possible early jet break around T * + 670 s in the light curve (LC) at a significance level of 9.2σ, where T * = T 0 + 226 s after the the Fermi Gamma-ray Burst Monitor trigger time T 0 at 2022 October 9 at 13:16:59 UT (Veres et al. 2022;LHAASO Collaboration et al. 2023).This compelling finding suggests that VHE radiation of GRB 221009A emanates from an exceedingly narrow jet.However, our previous work (Ren et al. 2023), along with a series of pertinent studies (Gill & Granot 2023;Kann et al. 2023;Sato et al. 2023;O'Connor et al. 2023), has revealed that the multiband afterglow cannot be satisfactorily explained by a simple narrow top-hat jet.Such an oversimplified model leads to untimely fading of afterglow flux and mismatches in decay indices.Consequently, the adoption of a structured jet model becomes indispensable to explicate the data.

Jet Structure
Numerical simulations of relativistic jet propagation through the progenitor star's envelope (e.g., Geng et al. 2019;Gottlieb et al. 2021Gottlieb et al. , 2022) ) have revealed a rich diversity of angular and axial structures in outflowing matter (e.g., Salafia et al. 2020).These structures are commonly categorized into three distinct components: the jet, the jet-cocoon, and the ejecta-cocoon.Often, structured jets with power-law or Gaussian angular profiles are adopted to model these components (e.g., Ren et al. 2020).Additionally, in modeling certain GRB afterglow observations, a commonly used approach is the so-called two-component jet model (e.g., Huang et al. 2004;Filgas et al. 2011).In this study, we describe the jet as a tophat inner 'core' encompassed by a power-law structured outer 'wing' in the surrounding region.We consider an axisymmetric structured jet as where E k,iso (θ) = 4πε(θ) is the isotropic kinetic energy of jet with ε(θ) being the energy per unit solid angle, Γ 0 (θ) is the initial bulk Lorentz factor and θ j is the angle of edge of the jet component, respectively.Subscripts '1' and '2' are used to differentiate between the core and wing.θ c,2 represents the generalized core opening angle of the wing, which describes the angle at which the properties of the wing undergo a transition.We artificially set a minimum value of 1.4 for Γ 0 (θ).Using the aforementioned formula, we can mimic various jet structures, thus inferring possible formation mechanisms for the jet structure of GRB 221009A.

Circum-Burst Environment
The circum-burst environment of GRBs exhibits a complex and diverse nature.Alongside the conventional environments of free stellar wind and interstellar medium (ISM), abrupt density transitions in the surrounding medium are often mentioned (Dai & Lu 2002).Some studies have explained the sudden re-brightening, fluctuations, and changes in decay slopes in the afterglow LCs of GRBs using models incorporating density jumps in the medium (e.g., Dai & Wu 2003;Monfardini et al. 2006;Li et al. 2020a).Numerical simulations have shown the mass outflows from massive stars could create a stellar wind bubble within a certain radius surrounding the star and eventually transition to a uniform ISM (e.g., Eldridge et al. 2006).
Initial researches did not reach consensus on the environment of GRB 221009A, e.g., Kann et al. (2023), O'Connor et al. (2023), andSato et al. (2023) preferred homogeneous, Laskar et al. (2023) and Ren et al. (2023) preferred wind-like.Subsequently, LHAASO reported the detected VHE afterglow with two distinct ascending slopes of ∝ t 15 and ∝ t1.8 , and they argue that the slope of ∝ t 15 could be the result of early energy injection.Furthermore, the observed slope of ∝ t 1.8 deviates from the typical wind environment but aligns with the scenario of a homogeneous medium surrounding the burst (LHAASO Collaboration et al. 2023), indicating a discrepancy to the free stellar wind environment implied by the later multiband afterglow observations (Gill & Granot 2023).In order to reconcile this problem, we propose an 'inverted' model that assumes a homogeneous environment near the progenitor star and transitions to a wind-like environment at a specific radius r tr , where A ⋆ = Ṁ /10 −5 M ⊙ yr −1 4πv w /1000 km s −1 −1 is the so-called wind parameter.

THE ACQUISITION AND USAGE OF DATA
We summarize here the database we used during this work.We use the VHE LC data and SED data reported in LHAASO Collaboration et al. (2023) 1 .We processed the publicly available data from Fermi-LAT and Swift-XRT and applied them in the fit.The optical data were obtained from Table 5 column SF11 (i.e., the Galactic extinction is E(B − V ) Gal = 1.322 ; Schlafly & Finkbeiner 2011) of Kann et al. (2023).We applied host extinction corrections to the data in the r and z bands.Due to uncertainties in the extinction factors, different articles provide varying values.We considered an appropriate extinction correction as E(B − V ) host = 0.3 with R V = 2.93 (Pei 1992) on this basis, but did not consider the subtraction of late supernova and host galaxy contributions (e.g., Blanchard et al. 2023;Fulton et al. 2023;Levan et al. 2023;Shrestha et al. 2023;Srinivasaragavan et al. 2023), so the late optical data may exceed the model curves.We took radio data from three articles, Bright et al. (2023), Laskar et al. (2023), andO'Connor et al. (2023).The radio LCs used the data from the first two papers (Figure 3), and radio spectra from all of them (Figure 5).

FITTING METHOD
In order to fit GRB 221009A, we used the ASGARD package to generate multiband LCs and spectral energy distributions (SEDs).Package details are presented in the Appendix A. In this work, we have numerically solved both forward shock (FS) and reverse shock (RS) dynamics (see A.2.1).We solve a series of evolution equations independently for the outflow at different jet angles, which includes the assumption that there is no causality between the elements at different angle, and leads to the result of structured forward and reverse shock radiation.The time-dependent electron continuity equation is solved with first-order accuracy which includes the Compton parameter that is associated with the energy of electrons.We accounted for synchrotron radiation and synchrotron self-Compton (SSC) emission, incorporating the corrections due to the Klein-Nishina effect.Additionally, we considered γγ annihilation effects, the equal-arrival-time surface effect, and the extragalactic background light (EBL) absorption.The EBL model we used is taken from Saldana-Lopez et al. (2021) We perform posterior parameter inference using the Markov Chain Monte Carlo (MCMC) technique, employing the Python package emcee as the sampler (Foreman-Mackey et al. 2013).The sampling was carried out in a 22 dimensional parameter space.The model parameters include where the subscription 'f ' and 'r' are used to differentiate between FS and RS, '1' and '2' mark the core and wing component of jet, respectively.
For simplicity, we have disregarded RS of the core component and focused solely on the radiation from the FS, as there is a lack of very early radio-to-optical band observations to say the RS exist or not.In contrast, we considered both FS and RS radiation in the wing of the jet, as early radio observations exhibited signs of RS emission (Bright et al. 2023).In our modeling, the electron acceleration efficiency ξ e,r,2 = 1 of RS is fixed due to the lack of sufficient data to constrain this parameter effectively.
In our fitting process, we simultaneously incorporated both the radio-to-TeV LCs data and the GeV-to-TeV SEDs data.This stringent combination provides a highly constraining posterior distribution for the model parameters.Our fitting procedure commences with a uniform prior distribution for the parameters.We adopted a uniform error σ of 6.03 +0.24   10% for all observed data points to ensure equal weighting of the data.Additionally, we did not account for potential systematic errors that might arise from differences between instruments3 .
We sampled the resulting posterior probability density using 50 walkers in emcee, which we ran over 7000 iterations.The resulting marginalized posterior probability density distributions, after discarding the initial 5000 iterations of the chain as burn-in, are shown in Figure 1.The obtained values of parameters are reported in Table 1, we have categorized the parameters into five groups based on their roles: core parameters, wing parameters, wing RS parameters, jetstructure parameters, and environment parameters.

Basic results
We find that the core component exhibits a typical set of afterglow parameters, with a large isotropic energy E k,iso,1 = 2.1 × 10 55 erg and a typical initial bulk Lorentz factor Γ 0,1 = 178.Thus the efficiency of prompt emission to the total energy is η γ = 42% (An et al. 2023;Yang et al. 2023).We calculate the beaming-corrected jet energy E k,j,1 ≈ 2.75 × 10 50 erg for the core, and total energy E k,j,tot ≈ 1.71 × 10 51 erg for whole jet.Fulton et al. (2023) has inferring an explosion energy of possible exist supernova SN 2022xiw as E k,SN ≈ 3−9×10 52 erg.Also, Srinivasaragavan et al. (2023) inferred another possible explosion energy E k,SN ≈ 1.6 − 5.2 × 10 52 erg.Thus, we find the allocation of burst energy is η = E k,j,tot /(E k,j,tot + E k,SN ) = 1 ∼ 10%, which consistent with the founding of Lü et al. (2018).We find that Γ 0,1 deviates slightly from the 2σ region of Γ 0 − E γ,iso correlation obtained by Liang et al. (2013), but is consistent with those found by the fittings from prompt emission of Yang et al. (2023).Compared to LHAASO Collaboration et al. (2023), we find that p f,1 = 2.345 is well consistent with the observed spectral index, but a smaller half-opening angle θ j,1 = 0.3 • is obtained.Although such a small half-opening angle presents challenges for the jet production and collimation process, we still observe a consistent correlation between the bulk Lorentz factor and the half-opening angle of jet θ j,1 ∼ 1/Γ 0,1 , satisfying the requirements of causality and collimation, indicating a selfconsistent parameter combination.We find the wing component operates in the mildly relativistic regime with bulk Lorentz factors in the range of tens, meanwhile the inferred half-opening angle is small, θ j,2 ∼ 1.3 • .We present the inferred angular-dependent isotropic energy E k,iso (θ) and bulk Lorentz factor Γ 0 (θ) distributions of the jet in Figure 2.

Jet Structure
The structure of GRB jets has been a topic of considerable interest.When studying the afterglow of the short GRB 170817A, the jet of GRB 170817A has been widely modeled as either Gaussian or power-law profiles, yielding exciting successes (e.g., Ren et al. 2020).A thought-provoking question arises as to whether the jet structure of long and short GRBs can be consistently described using the same form (Salafia et al. 2020).Although some numerical simulations have shown commonalities between them, observational constraints on the angular-dependent structure of jet has long been elusive.
Figure 2 shows the distribution of isotropic kinetic energy and bulk Lorentz factor with angle of GRB 221009A.The findings from E k,iso (θ) are highly intriguing.Despite being results obtained through free parameters fitting, it reveals a remarkably smooth transition between the core and wing regions.Unlike the energy, the bulk Lorentz factor exhibits a steep transition.The core region exhibits lower baryon loading, resulting in faster acceleration and dominating the emission of TeV afterglows (Figure 3).The wing component operates in the mildly relativistic regime Γ 0 ∼ 50 and gradually transitions to Γ 0 ∼ 20 at larger angles.This finding suggests a fundamentally different physical origin between the core and wing components.As shown in Figure 2, the overall angular structure distribution of the jet is inconsistent with a two-component jet.Instead, this structure is more reminiscent of the properties exhibited by jetcocoon systems in numerical simulations (Gottlieb et al. 2021(Gottlieb et al. , 2022)).Another evidence is the small wing half-opening angle θ j,2 ∼ 1.3 • .Therefore, in terms of definition, the wing component is closer to a jet-cocoon.We notice that our study lacks modeling of the ejecta-cocoon with lower velocity, which may be the reason for inability to explain the low-frequency radio observations satisfactorily (Laskar et al. 2023).Gottlieb et al. (2021) proposes the presence of a jet-cocoon interface region in the angular structure of GRB jets, where intense mixing processes result in a smooth transition of jet properties from the core to the wing regions.We observe similar energy characteristics of GRB 221009A, but the sharp change in the bulk Lorentz factor reveals a E n e r g y ( e V ) lack of effective material mixing mechanisms in the interface region, which contrasts significantly with hydrodynamic simulations.As a result, our findings may suggest the magnetized nature of the jet, or, reflect the subsequent impact of the continuous jet from central engine on the jet structure after the jet-head breaking out of the envelope.The structure of jet constrained from GRB 221009A differs from that of GRB 170817A.First, the jet of GRB 170817A has a core with half-opening angle θ c ∼ 3 • − 5 • and a broad wing θ j ∼ 20 • (e.g., Ren et al. 2020).GRB 221009A obviously has smaller values.Then, the index of decayed energy of GRB 221009A in the wing region E(θ) ∝ θ −5.3 appear to be consistent with those of GRB 170817A (e.g., E(θ) ∝ θ −5.5 , Ghirlanda et al. 2019; E(θ) ∝ θ −4.5 , Hotokezaka et al. 2019; E(θ) ∝ θ −6 , Ren et al. 2020).But, although there is a jump in the bulk Lorentz factor between the two regions of the jet of GRB 221009A, the trend of bulk Lorentz factor Γ(θ) ∝ θ −2 seems more flat than the jet of GRB 170817A (e.g., Γ(θ) ∝ θ −3.5 , Ghirlanda et al. 2019 Ren et al. 2020).It may imply a common physical origin of the structures of both jets, i.e., the jet passing through a dense material envelope, but there is a significant difference in the nature of the envelope that the two jets encounter.GRB 221009A, being a long GRB, has its jet passing through a dense and static stellar atmosphere.In contrast, the jet of short GRB 170817A encounters an expanding and relatively thin ejecta during the propagation.As a result, the jet of GRB 221009A undergoes a stronger interaction with the surrounding material and experiences highly collimation.Besides, the mixing of material in the cocoon maybe more efficiency, so that the Lorentz factor more evenly distributed but has lower values.We believe this explanation could providing a reasonable interpretation for the unique jet structure observed in GRB 221009A.
Our findings also suggest that a subclass of low-luminosity GRBs with TeV observations may originate from such jet structures.The axis of jet deviates slightly from the line-of-sight, and the observed afterglow could arise from the radiation of a slow-moving jet-cocoon, as seen in events such as GRBs 190829A and 201015A (Zhang et al. 2021;Salafia et al. 2022;Zhang et al. 2023b).

5.3.
Comparing ϵ e , ϵ B , and ξ e of Core and Wing Components For the core component, the derived values of ϵ e,f,1 = 6.6 × 10 −2 and ϵ B,f,1 = 9.3 × 10 −7 are in good agreement with typical values of GRB samples (Barniol Duran 2014; Wang et al. 2015;Beniamini & van der Horst 2017;Duncan et al. 2023).The small value of ϵ B,f,1 suggests that turbulence behind the shock amplifies the magnetic field with very low efficiency, or the turbulence rapidly decays with increasing distance from the shock front (Lemoine 2013(Lemoine , 2015)).We F r e q u e n c y ( H z ) note that the electron acceleration efficiency parameter, ξ e,f,1 = 0.44, indicates a considerable fraction of electrons are efficiently accelerated to a relativistic distribution.This requires the core component of the jet to possess a powerful acceleration mechanism of particles.
The parameters for the wing component are relatively more intriguing.Compared to the core, the wing component has a smaller ϵ e,f,2 = 4.6 × 10 −3 , larger ϵ B,f,2 = 1.1 × 10 −3 , and ξ e,f,2 = 7.6 × 10 −3 that is lower by almost two orders of magnitude.Considering that both components exist in the same surrounding medium and have similar isotropic kinetic energy of jet, we can speculate that the lower electron acceleration efficiency is associated with smaller bulk Lorentz factor and a larger ϵ B .Interestingly, TeV-afterglow-detected nearby low-luminosity GRB 190829A (H.E. S. S. Collaboration et al. 2021) seems to provide evidence for this relevance.GRB 190829A has a moderate bulk Lorentz factor Γ 0 ∼ 45 − 55 (Zhang et al. 2021;Salafia et al. 2022;Huang et al. 2023).As demonstrated by Salafia et al. (2022), it exhibits a very low electron acceleration efficiency ξ e,f = 7 × 10 −3 (wide prior).Although there are differences in the fitting results for other parameters, this relevance may indicate a direction for future research.Further investigations with a larger sample of well-fitted cases are expected to verify this correlation.

Fitting Results of LCs and SEDs
In Figure 3, we present the multiband afterglow LCs obtained using the best-fit parameters, with different line styles distinguishing the emissions from core FS, wing FS, and wing RS, respectively.Furthermore, in Figure 4, we display the fitted results for the GeV-to-TeV energy range SEDs, with different line styles representing the core and wing components.Lastly, Figure 5 demonstrates the fitted ν − F ν spectra for the radio-to-X-ray bands.
The TeV emission is found to be dominated by the core component, as demonstrated in Figure 3.By assuming our medium model (Equation 3), all LHAASO observations can be consistently explained.The rising phase of LC is caused by the uniform ambient medium environment at r < r tr .The rapid decline of LC, observed around T * + 700 s, is attributed to a smooth jet break occurring within the free stellar wind environment.
The broadband emission, spanning from radio to X-ray, is primarily determined by the radiation from both FS and RS of the wing component.Observations in the optical and X-ray bands are dominated by FS radiation.Some detectors indicate the presence of persistent power-law components in the spectra of tailed prompt emission.The inferred X-ray LC is extend from T * + 100 to T * + 2000 s, showing a transition from flat to steep decay.This radiation component is attributed to external shock emission in some studies (An et al. 2023;Lesage et al. 2023;Zhang et al. 2023a).However, we suggest that this early X-ray LC cannot be explained by external shock radiation because a satisfactory simultaneous fit cannot be reached for both the radio and X-ray bands.Instead, we are inclined to attribute this radiation to the 'late prompt emission'.In other words, the continuous and smooth emission characterized by a power-law decay over time is produced by the internal dissipation of shells from late-time central engine activity.Meanwhile, the transition in the decay slope can also be explained self-consistently within this model (Ghisellini et al. 2007).
The early-time radio observations show a good fit with both FS and RS radiation, whereas RS radiation dominates the late-time radio observations.However, we still lack a precise explanation for the late-time radio observations.Specifically, the low-frequency flux is significantly underestimated (see Figure 5).It is possible that we have overlooked an additional trans-relativistic outflow component known as the ejecta-cocoon, given our employed modeling strategy.Introducing this component into the model could potentially lead to an improvement in the fitting.However, considering the satisfactory results achieved with our current model, which already includes up to 22 parameters, it may not be worthwhile to add several additional parameters to describe the properties of the ejecta-cocoon.It would be more advantageous to leave this aspect of modeling for future investigations, as continuous radio observations will provide late-time data.Hence, we only briefly discuss this possibility here.
Another potential explanation for the late-time radio emission is the reheating of low-energy electrons through the synchrotron self-absorption (SSA) process.This results in the accumulation of heated electrons around the selfabsorption frequency, leading to a noticeable bump in the SEDs near this frequency, as demonstrated in previous studies (e.g., Yang et al. 2016;Li et al. 2020b).Additionally, the investigation of thermal electron radiation presents an interesting direction, as it can also have an impact on the radio band (Eichler & Waxman 2005;Fukushima et al. 2017;Margalit & Quataert 2021;Warren et al. 2022;Nedora et al. 2023).However, both SSA and thermal electron radiation are only significant during the early stages of afterglow.Consequently, their effectiveness in explaining the late-time low-frequency radio emission might be limited.
During the investigation, we also observed that the behavior of GRB 160625B at low frequencies during the afterglow phase is similar to that of GRB 221009A (Alexander et al. 2017;Cunningham et al. 2020;Kangas et al. 2020).Notably, at 12 and 22 days, GRB 160625B exhibited a re-emergence of low-frequency radio emission.This phenomenon cannot be explained by the RS mechanism and is consequently attributed to an extreme scattering event (ESE) (Alexander et al. 2017).In contrast to GRB 160625B, GRB 221009A has a low Galactic latitude.This factor increases the probability of encountering dense plasma structures that function as lenses along the path of photon propagation (Tiengo et al. 2023;Vasilopoulos et al. 2023).Consequently, GRB 221009A might be another example of an ESE.

Circum-Burst Environment
The environment surrounding long GRBs has garnered considerable research interest (Mirabal et al. 2003;Chevalier et al. 2004;Ramirez-Ruiz et al. 2005;Eldridge et al. 2006).Although observations have confirmed the association of long GRBs with supernovae, indicating their connection to the collapse of massive stars (e.g., Woosley & Bloom 2006), studies of long GRB afterglow samples reveal that a significant portion of afterglows does not exhibit characteristics of a free stellar wind environment (e.g., Schulze et al. 2011;Liang et al. 2013;Yi et al. 2013Yi et al. , 2020;;Tian et al. 2022).This complexity suggests the diverse origins of the environments surrounding the bursts.Intriguingly, evidence of one or multiple dense material shells surrounding the progenitor star has been discovered in certain supernovae (e.g., Corsi et al. 2014;Milisavljevic et al. 2015) and interacting superluminous supernovae (e.g., Liu et al. 2018).The density distribution of these shells often suggests a homogeneous medium, as opposed to a free stellar wind medium.In summary, these findings imply complex and unstable material ejection behavior during the post-main-sequence lifetime of massive stars (Langer 2012).
This work utilizes a phenomenological model described by Equation 3 to tackle the challenge of explaining the VHE afterglow observations in GRB 221009A.By considering the transition from a homogeneous to a wind-like medium at r tr = 3.7 × 10 16 cm, we have achieved a successful fit to the data.The number density of medium within r ⩽ r tr is n tr = 1125 cm −3 .It should be noted that the VHE LC exhibits a rapid rise at approximately T b ∼ T * + 5 s.In our proposed scenario, the rapid rise can be attributed to the initiation of an external shock when the jet collides with the shell located at r s .This collision implies a lower density at r < r s compared to n tr .Based on this assumption, we can estimate the position of the shock front at this time to be r s = 2Γ 2 0,1 cT b /(1 + z) = 8.2 × 10 15 cm.In this context, some TeV photons detected prior to T * + 5 s, surpassing the background level, can be attributed to either internal shock events or collisions between the jet and material clumps or shells located at r < r s (Marchenko et al. 2007;Cantiello et al. 2009).The lower energy gamma-rays may be obscured by the prompt emission during this period.
Our modeling of the VHE afterglow of GRB 221009A suggests the presence of a homogeneous shell between r s and r tr .Assuming a spherically symmetric distribution of the shell, we calculate its mass to be approximately 2×10 −4 M ⊙ , which is smaller compared to the shell masses observed in supernovae and superluminous supernovae.In addition, if we assume a stable free stellar wind, this shell cannot be attributed to the jet shell of the precursor of GRB 221009A.If it were, it would gather considerably more material than our calculations show, contradicting our results.
We propose that this shell reflects the evolution of material outflow just prior to the death of the star.A sudden decline in the strength of the stellar wind, potentially accompanied by a minor ejection event, may explain the presence of this shell.Taking into account the assumption that the progenitor is a typical Wolf-Rayet star, the velocity of the ejection shell is often slower than the stellar wind velocity (Milisavljevic et al. 2015).Consequently, the structure of the ejection shell becomes dispersed during its propagation, thus preventing the formation of shocks in the free stellar wind.This is a plausible explanation for the smooth transition in the density of the medium.

SUMMARY AND CONCLUSION
GRB 221009A, as a powerful GRB event, has provided a rich and complete set of GRB afterglow observations, offering a comprehensive showing stage for the decades-developed theories in the field of GRB afterglow.Our study showcases the success of theoretical efforts, but also show some outstanding questions in future researches.
• Our study suggests that GRB 221009A possesses a core+wing structure, where E k,iso (θ) shows a smooth transition, but Γ 0 (θ) displays a clear interface stratification between the core and wing.We neglected the ejecta-cocoon component in our modeling, but the underestimation of late-time low-frequency radio flux suggests its possible existence, and future observations of GRB 221009A will help answer this question (Bright et al. 2023;Laskar et al. 2023).The structure of the jet in GRB 221009A is differ from that of short GRB 170817A, possibly due to differences in the jet magnetization and the envelope properties of progenitor.Our findings also suggest that a subclass of low-luminosity GRBs with VHE observations, e.g., GRBs 190829A and 201015A, may originate from such jet structures.
• The fittings demonstrates that SSC emission from core component can successfully explain the LHAASO observations, while the lower energy afterglow requires the presence of the synchrotron emission from wing component for explanation.However, the previously reported > 10 TeV photons clearly cannot be explained by the SSC process (Huang et al. 2022).Further investigation is needed to determine whether observable hadronic processes are present during the afterglow phase (e.g., Sahu et al. 2022;Isravel et al. 2023).
• Our fitting reveal a distinct transition in the burst environment from uniform to wind-like medium, indicating the presence of intricate pre-outburst mass ejection processes.The detection by LHAASO of the early VHE afterglow provides a novel avenue for investigating the burst environment of progenitor star.More observations of VHE-afterglow-detected GRBs will help to unveil more details about the burst environment.
• By comparing the FS parameters of the core and wing components, we suggest a possible correlation between electron acceleration efficiency ξ e and both the Lorentz factor of the shock and the magnetic field equipartition factor ϵ B .To test this hypothesis, analyses require using a set of well-sampled afterglow data of GRBs, and performing plasma numerical simulations under extreme conditions.
• The complex radio and millimeter emission are difficult to explain within the scope of standard synchrotron emission theory (Bright et al. 2023;Laskar et al. 2023;O'Connor et al. 2023).Although we did not discuss the aspect of evolving microphysical parameters, ϵ e , ϵ B , ξ e , and p, in our work, considering such evolution may potentially improve the current fitting results (e.g., Yang et al. 2018).

APPENDIX
A. MODEL DETAILS This section describes the modeling details of the code package utilized for fitting.

A.1. Code Framework of ASGARD
A Standard Gamma-ray Burst Afterglow Radiation Diagnoser (ASGARD) is a software package designed to generate multiband afterglow lightcurves (LCs) and spectral energy distributions (SEDs) for specific events.The package provides a Python interface with underlying FORTRAN code that allows users to easily conduct research after a simple learning and setup process.
To ensure code reusability, extensibility, and readability, ASGARD follows a modular approach where different modules handle specific calculations.Python code is used to call these modules.ASGARD currently encompasses three main components: jet dynamics, electron distributions, and radiation physics.Radiation physics includes synchrotron radiation and synchrotron self-Compton (SSC) radiation of electrons.
All the modeling in ASGARD is based on detailed numerical calculations, aiming to realistically describe the radiation behavior of the afterglow.While ASGARD follows the 'standard' fireball model, it also incorporates post-standard models developed over decades, such as energy injection and structured jets.The extensible architecture of the package allows for the inclusion of additional considerations in the future.A brief description of the numerical methods and models used is provided in the subsequent subsections.

A.2. Dynamics of Jet propagation
During the expansion of a relativistic fireball in a surrounding medium, two shock waves are created: a forward shock (FS) and a reverse shock (RS).Upon the formation of shocks, a region of shock interaction arises, characterized by two distinct components: the shocked medium and the shocked shell.These two regions are separated by the contact discontinuity surface.The evolution of these shock waves and the physical properties of the fireball involves solving a set of equations that couple relativistic (magneto)hydrodynamic waves.Blandford & McKee (1976) derived a self-similar solution for adiabatic relativistic hydrodynamic waves, which provides insights into the time-dependent evolution of the fireball shocks.However, to obtain a more precise description of the evolution of jet, a set of differential equations incorporating radiation cooling in one dimension is employed.These equations account for the energy lost through radiation as the fireball interacts with its surroundings.
The solution of these equations relies on the assumption of the single-shell hypothesis where the majority of the mass and energy of the fireball are concentrated within a small radius in the vicinity of the shock region.This assumption simplifies the calculations and allows for a focused analysis of the shock dynamics and radiation behavior.

A.2.1. Dynamics of coupled forward-reverse shock
At the beginning of external shock generation, RS propagates into the jet shell, and FS/RS dynamics are coupled at this stage.There is a possibility that the RS strength may be suppressed if the jet is magnetized (e.g., Zhang & Kobayashi 2005;Chen & Liu 2021).At present, ASGARD does not include magnetizing jet modeling.As a substitute, we used model of Yan et al. (2007) to analyse RS dynamics in a pure fireball scenario.
A homogeneous cold shell is expelled from the central engine, carrying an isotropic energy E k,iso , and possesses a width ∆ 0 .As it expands relativistically with a Lorentz factor Γ 0 , the shell mass is given by m 0 = E k,iso /Γ 0 c 2 , where c is the speed of light.
Within the shocked region, it is assumed that the bulk velocity and energy density are homogeneous, implying that γ 2 = γ 3 and e 2 = e 3 .Hereafter, the subscript '2' designates the shocked medium, the subscript '3' corresponds to the shocked jet shell, and the subscript '4' pertains to the unshocked jet shell, i.e., γ 4 = Γ 0 .

A.3. Solving of the Electron Continuity Equation
We denote the instantaneous electron distribution as dN e /dγ ′ e , of which the evolution can be solved based on the continuity equation of electrons (Fan et al. 2008) where "′" marks the co-moving frame of shock, and  (Sari et al. 1998), with Γ = γ 2 for FS, and Γ = γ 34 for RS, respectively.We assume that a fraction ξ e of the electron population is shock-accelerated (Eichler & Waxman 2005).Evidently, the overall inflow of nonthermal electrons should be The maximum Lorentz factor of electron is γ ′ M = 9m 2 e c 4 /[8B ′ q e 3 (1 + Y )] with q e being the electron charge (Kumar et al. 2012), where Y is the Compton parameter.We note that the Compton parameter has been solved based on the work of Fan & Piran (2006) in their appendix, ϵ e and ϵ B are the equipartition factors for the energy in electrons and magnetic field in the shock, respectively.For the photons with frequency higher than ν′ , the Compton parameter should be suppressed significantly since it is the Klein-Nishina regime, where ν′ is governed by γ ′ e hν ′ ∼ Γm e c 2 .Since we have consider the case of p > 2, for the slow-cooling case Ren et al.
A12) also, for the fast-cooling case (γ ′ m > γ ′ c ), The electron populations in both FS and RS has treated using the method described above.For electrons in RS, n(r) in Equation A10 should be considered as n(r) = n 4 .Unlike the continuous injection of electrons in the case of FS, no new electrons are injected after the crossing time t ′ × of RS which has traversed all the matter in the jet shell, completely replacing zone 4 with zone 3.

A.4. Radiation Mechanism
In the X-ray/optical/radio bands, the main radiation mechanism of the electrons in GRB jets is synchrotron radiation (Sari et al. 1998;Sari & Piran 1999).The emitted spectral power of synchrotron radiation at a given frequency ν ′ of a single electron is where we have assumed that the direction of the magnetic field is perpendicular to the electron velocity.
is the modified Bessel function of 5/3 order, and ν ′ c = 3q e B ′ γ ′ e 2 /(4πm e c), respectively.Thus, the spectral power of synchrotron radiation of electrons dN e /dγ ′ e at a given frequency ν ′ is The photon seed spectra of synchrotron radiation is then be calculated as The emission of the SSC process is calculated based on the electron spectrum dN e /dγ ′ e and target seed photons of synchrotron radiation n ′ γ,syn (ν ′ t ) from the synchrotron radiation (e.g., Fan et al. 2008;Nakar et al. 2009;Geng et al. 2018) e hν ′ t /m e c 2 , and w = hν ′ /γ ′ e m e c 2 , respectively.One can also derive the photon seed spectra of SSC radiation The total spectral power and photon seed spectra then be P ′ tot (ν ′ ) = P ′ syn (ν ′ ) + P ′ SSC (ν ′ ) and n ′ γ,tot = n ′ γ,syn + n ′ γ,SSC , respectively.We also considered the γγ annihilation effects (e.g., Gould & Schréder 1967;Geng et al. 2018;Huang 2022).Since the cross section of γγ annihilation reads ) and where β cm is the dimensionless velocity of electron (positron) in the center-of-momentum frame, µ = cos θ with θ being the angle of collision photons, and ν ′ t is the frequency of target photons from synchrotron and IC radiation.It is obvious that there is a physical value only when hν ′ hν ′ t (1 − µ) ⩾ 2 m e c 2 2 .Thus, the optical depth for a gamma-ray photon with frequency ν ′ is expressed as (e.g., Murase et al. 2011) The synchrotron self-absorption (SSA) effect is also considered in our numerical calculations.The optical depth is the function of both the electron distribution and synchrotron radiation power of single electron, note that γγ annihilation and SSA dominate the absorption of high-energy and low-energy photons, respectively.
A.5.The Geometric and Observational Effects GRBs are believed to be generated by the motion of ultra-relativistic jets within a half-opening angle θ j .Due to geometric and relativistic Doppler effects, the observed radiation spectrum cannot be simply reproduced by the intrinsic spectrum, i.e., P ′in tot (ν ′ ).In this work, we set the GRB jet as an on-axis-observed jet, i.e., θ v = 0. Assuming that the observed frequency is ν obs , the frequency transformed to the comoving frame is denoted as ν ′ (ν obs ) = (1 + z)ν obs /D, where z is the redshift, (A24) where D L is the luminosity distance from burst to the Earth.The integration of θ is performed over an elliptical surface that photons emitted from the surface have the same arrival time for an observer (Geng et al. 2018), Gamma-ray photons traveling through the universe may undergo γγ annihilation with extragalactic background light (EBL) photons, resulting in additional absorption of high-energy photons.When considering the EBL absorption, the finally observed flux reads F obs tot (ν obs ) = F in tot (ν obs ) exp −τ EBL (ν obs , z) .(A26) A.6.Numerical Method The evolution of dynamics is a system of ordinary differential equations which is solved using the fourth-order Runge-Kutta method.For the distribution equation of electrons, which is a partial differential equation, a fully implicit first order scheme is implemented (Chang & Cooper 1970;Chiaberge & Ghisellini 1999;Huang 2022), with grid transformations to improve computational accuracy (e.g., Geng et al. 2018).During computations related to the radiation spectrum, the Euler method and linear interpolation are used.Therefore, the computed results used for fitting are believed to have first-order accuracy.
In this study, we assume that the jet is structured.This implies that the development of the jet at different angles is not connected.Consequently, we calculate the intrinsic radiation spectrum for each angle independently, and the final result is obtained by using the interpolation technique.We do not take into account any potential lateral expansion effects.This simplification is suitable when the jet is highly relativistic; however, it may lead to an overestimation of the flux after the jet break time.

Figure 1 .
Figure 1.Corner plot of the MCMC posterior sample density distributions.

Figure 3 .
Figure 3. Multiband afterglow light curves plotted with best fitting result.Filled circles are data points and open triangles are upper limits.Bands are in different colors with shift factors in legend.Solid lines represent the total emission, dashed lines represent core forward shock emission, dotted lines represent wing forward shock emission, and dash-dotted lines represent wing reverse shock emission, respectively.

Figure 4 .
Figure 4. SEDs plotted with best fitting result for GeV-to-TeV observations, data were obtained from LHAASO Collaboration et al. (2023) and Ren et al. (2023).Time slices are in different colors with shift factors in legend.Solid lines represent the total emission, dashed lines represent the synchrotron emission, and dotted lines represent the SSC emission, respectively.

Figure 5 .
Figure5.SEDs plotted with best fitting result for radio-to-X-ray observations.Time slices are in different colors with shift factors in legend.To distinguish between the datasets, the dashed lines represent the data fromLaskar et al. (2023), while the solid lines correspond to the data from O'Connor et al. (2023).We interpolated data fromBright et al. (2023) to convert it to 0.25 day.
A9) is the cooling term which has included the cooling effects of synchrotron radiation, SSC process, and adiabatic expansion.The swept-in electrons by the shock are accelerated to a power-law distribution of Lorentz factor γ ′ e , i.e., Q ∝ γ ′ e −p for γ ′ m ⩽ γ ′ e ⩽ γ ′ M , where p(> 2) is the power-law index.Here, γ ′ m = ϵ e /ξ e (p − 2)m p Γ/[(p − 1)m e ] + 1 A13) where γ′ e = 4πm e cν ′ /3q e B ′ , and γ ′ c = 6πm e c/(σ T ΓB ′ 2 t ′ ) is the efficient cooling Lorentz factor of electrons with σ T being the Thomson scattering cross section.
intrinsic spectral power of afterglow is β cos Θ) is the Doppler factor, and Θ represents the angle between the direction of motion of a fluid element in the jet and the line-of-sight.If we take the equal-arrival-time surface (EATS) effect into account(Waxman 1997), the observed intrinsic spectral flux should beF in tot (ν obs ) tot [ν ′ (ν obs )] D 3 sin θ 2 dθ,

The
Rainbow Bridge of Asgard in Norse mythology, has inspired R.J. to name the ASGARD package.It traverses the universe with shining, just like the GRB afterglow.This work was supported by the National Natural Science Foundation of China (grant No. 11833003).This work used data and software provided by the Fermi Science Support Center and data supplied by the UK Swift Science Data Centre at the University of Leicester.Software: Matplotlib (Hunter 2007), NumPy (Harris et al. 2020), emcee (Foreman-Mackey et al. 2013), corner (Foreman-Mackey 2016), Astropy (Astropy Collaboration et al. 2013), extinction (Barbary 2016), SciPy (Virtanen et al. 2020), GetDist (Lewis 2019) , same as what LHAASO Collaboration et al. (2023) used.
is the dynamical timescale of syn-