MAGIC detection of GRB 201216C at $z=1.1$

Gamma-ray bursts (GRBs) are explosive transient events occurring at cosmological distances, releasing a large amount of energy as electromagnetic radiation over several energy bands. We report the detection of the long GRB~201216C by the MAGIC telescopes. The source is located at $z=1.1$ and thus it is the farthest one detected at very high energies. The emission above \SI{70}{\GeV} of GRB~201216C is modelled together with multi-wavelength data within a synchrotron and synchrotron-self Compton (SSC) scenario. We find that SSC can explain the broadband data well from the optical to the very-high-energy band. For the late-time radio data, a different component is needed to account for the observed emission. Differently from previous GRBs detected in the very-high-energy range, the model for GRB~201216C strongly favors a wind-like medium. The model parameters have values similar to those found in past studies of the afterglows of GRBs detected up to GeV energies.


INTRODUCTION
Gamma-ray bursts (GRBs) are sources exhibiting bright electromagnetic emission in two phases called prompt and afterglow.The former peaks at hard X-ray and soft gamma-ray energies, lasting between a fraction of a second and hundreds of seconds.In particular, the prompt temporal behavior shows short time scale variability down to milliseconds.Although its origin is not completely understood (for a review, see Kumar & Zhang 2015), recent evidence is pointing to a synchrotron origin (Zhao et al. 2014;Zhang et al. 2016;Oganesyan et al. 2017Oganesyan et al. , 2018Oganesyan et al. , 2019)).The afterglow radiation partly overlaps with the prompt and evolves over longer timescales, up to several months after the GRB onset.The emission in this phase decays smoothly with time as a power law and it can be detected in several energy bands, from radio up to gamma rays, and is interpreted as synchrotron and Inverse Compton emission mostly from electrons accelerated in the external shock (Sari et al. 1998;Panaitescu & Kumar 2000).
GRBs are classified as short and long depending on whether their duration in terms of  90 , the time interval containing 90% of the total photon counts, is shorter or longer than two seconds.While this observational definition is widely adopted, a more physical classification comes from the progenitor system at the origin of the bursts.In this context, short GRBs are thought to be produced as the result of the merger of binary systems of compact objects involving at least one neutron star (NS).The only confirmation of such association is the short GRB 170817A, which was detected in coincidence with a gravitational wave signal generated by a NS-NS merger (Abbott et al. 2017;Goldstein et al. 2017).On the other hand, long GRBs are often associated with supernovae of type Ib/c, when detectable (e.g. if redshift is  ≲ 1).The supernova emission peaks several days after the GRB onset, when it outshines the decaying optical afterglow of the burst itself (Woosley & Bloom 2006).
The afterglow phase of GRBs has been studied in detail over several wavelength bands thanks to numerous instruments both groundbased (covering the radio and optical wavelengths, and VHE gamma rays) and space-based (detecting X-rays and gamma rays).Such observations have made it possible to trace the origin of the multiwavelength afterglow emission to the synchrotron process (Mészáros 2002;Piran 2004).Such radiation is mostly produced by electrons accelerated at the so-called forward shock, when the GRB jet decelerates by interacting with the interstellar or circumstellar medium.Until recently, the afterglow was detected up to GeV energies by the Fermi-LAT instrument, with some hints of a possible tail extending to higher energies (Ackermann et al. 2014), where imaging atmospheric Cherenkov telescopes (IACTs) are more sensitive.The presence of emission in the very-high-energy (VHE,  > 100 GeV) range in the afterglow phase of GRBs was predicted, even before the operation of Fermi-LAT and IACTs, in several theoretical models involving either leptonic or hadronic processes.A breakthrough was achieved in 2019, when the detection of VHE emission in the afterglow of three long GRBs was reported.The MAGIC collaboration first reported the detection of GRB 190114C (Mirzoyan et al. 2019;MAGIC Collaboration et al. 2019a,b), followed by GRB 180720B and GRB 190829A detected by the H.E.S.S. telescopes (Abdalla et al. 2019; H. E. S. S. Collaboration et al. 2021).The detection of such sources with IACTs confirmed the presence of an emission in the VHE range.In particular, the spectral and temporal analysis of GRB 190114C showed that such emission is associated with a component, separate from the synchrotron one, well explained by synchrotron-self Compton (SSC) radiation from electrons accelerated at the forward shock.A similar conclusion can be drawn for GRB 180720B (see e.g.Wang et al. 2019), even though the multi-wavelength data available were not enough to perform a proper modeling.An unusual and controversial interpretation was put forward in the case of GRB 190829A.In H. E. S. S. Collaboration et al. (2021) the authors suggested that the emission could be attributed to a single synchrotron component, which extends over nine orders of magnitude in energy up to the TeV domain.This requires an acceleration mechanism that is able to overcome the limit resulting in the so-called burnoff limit for the energy of synchrotron photons (de Jager et al. 1996;Piran & Nakar 2010).
The studies on this small sample of events shows how the understanding of the afterglow phase in the VHE range is far from complete.Currently only a few events have a detection at VHE (or evidence, as in GRB 160821B, see Acciari et al. 2021), and different interpretations were proposed.However, the SSC scenario proved to be flexible and applicable to all the three GRBs detected at VHE.In order to investigate if such an interpretation may be universal to explain VHE afterglows, we present here the detection of the long GRB 201216C with the MAGIC telescopes.We use the available multi-wavelength data to model the broadband emission in the SSC scenario.We find that the SSC model provides a satisfactory interpretation of the MAGIC light curve and spectrum.
The paper is organized as follows.In Section 2, we summarise all the observations available for GRB 201216C.In Section 3 we discuss the MAGIC observations and data analysis.The results are presented in Section 4. In Section 5 we present the analysis of optical observations taken with the Liverpool Telescope and the other multiwavelength observations that we use to model the emission with a synchrotron and SSC scenario (discussed in Section 6).Finally, in Section 7 we summarise and discuss our findings.

OBSERVATIONS OF GRB 201216C
GRB 201216C was detected by Swift-BAT on December 16th 2020 at 23:07:31 UT (Beardmore et al. 2020) 1 , hereafter  0 .The burst was also detected by other space-based instruments including Fermi-GBM, ASTROSAT and Konus-Wind.The light curve by Swift-BAT shows a multi-peaked structure 2 from  0 − 16 s to  0 + 64 s, with a main peak occurring at ∼  0 + 20 s.
Observations at different times by the VLT, FRAM-ORM and the Liverpool Telescope confirmed the presence of the optical afterglow.The position of the optical counterpart is consistent with the refined position provided by Swift-XRT.VLT X-Shooter spectroscopy at ∼  0 +2.4 hours, covering the wavelength range 3200-22 000 Å, allowed the measurement of the redshift, estimated 3 to be  = 1.1.Based on the VLT photometry, the steep photon index of optical data suggests a significant extinction, making GRB 201216C a dark GRB (Vielfaure et al. 2020).
The afterglow was also detected in the X-ray band by Swift-XRT.The X-ray afterglow decay 4 can be described as a power law with temporal index  = 1.75 ± 0.09.
At higher energies, the burst was observed by HAWC starting at  0 + 100 s up to  0 + 3600 s, resulting in a non significant detection (Ayala 2020).
Detection of radio emission was reported by Rhodes et al. ( 2022) from 5 to 56 days after the burst, from 1 to 10 GHz.The radio flux at the time of detection is already decaying, although at a slow rate, except for the flux at 1 GHz, for which the flux is increasing between 30 and 40 days.
Finally, the burst was observed by the MAGIC telescopes in the VHE range.Details of such observations are given in the following section.

MAGIC OBSERVATION AND DATA ANALYSIS
MAGIC is a stereoscopic system of two 17-m diameter IACTs situated at the Observatory Roque de los Muchachos (ORM), La Palma, Canary Islands.For short observations, as the ones usually performed for GRBs, the integral sensitivity achieved by MAGIC in 20 min is about 20% of the Crab Nebula flux above 105 GeV for low zenith angles (see Aleksić et al. 2016 for details on the telescopes performance).
MAGIC received the alert for GRB 201216C at 23:07:51 UT ( 0 + 20 s) from the Swift-BAT instrument.The MAGIC telescopes automatically reacted to the alert and, after a fast movement, they reached the target at 23:08:27 UT ( 0 + 56 s).The observation was carried out in the so-called wobble mode around the coordinates provided by Swift-BAT, RA:01h05m26s Dec:+16d32m12s (J2000).In local coordinates, the observation started at zenith 37.1 • , lasting up to 01:30:08 UT reaching zenith 68.3 • .The weather conditions were very good and stable during all the data taking with a median atmospheric transmission value at 9 km a.g.l from LIDAR measurements of 0.96, with 1 being the transmission of a clear atmosphere (see Fruck et al. 2022;Schmuckermaier et al. 2023 for a description of the LIDAR instrument and correction of VHE data).The observation was performed under dark conditions.
MAGIC continued the observation on the second night for 4.1 h from  0 + 73.8 ks.The observational conditions were optimal with an average transmission above 0.9 at 9 km and dark conditions.The zenith angle changed from 17.0 • to 46.3 • with culmination at 11.7 • .The data on the second night was taken with the analog trigger system Sum-Trigger-II (described in Dazzi et al. 2021), which was not available during the first night of data taking.Sum-Trigger-II improves the sensitivity of MAGIC in the low-energy range below ∼ 100 GeV.In particular, the trigger efficiency, compared to the standard digital trigger, is two times larger for Sum-Trigger-II at 40 GeV.
The data analysis is performed using the standard MAGIC Reconstruction Software (MARS; Zanin et al. 2013).In order to retain as many low energy events as possible, an algorithm (Shayduk 2013;MAGIC Collaboration et al. 2020) where the calibration and the image cleaning are performed in an iterative procedure was adopted.This image cleaning was applied to the GRB data, gamma-ray Monte Carlo data, and to a data sample taken on sky regions without any gamma-ray emission (used for the training of the particle identification algorithm).Data analysis beyond this level is performed following the prescriptions described in Aleksić et al. (2016).The usage of Sum-Trigger-II, combined with the optimized cleaning algorithm, allows for a collection area an order of magnitude larger around 20 GeV when compared with the one obtained with the standard digital trigger.

RESULTS FROM THE VERY-HIGH-ENERGY DATA
In this section we show the results of the analysis performed on the data collected by MAGIC on GRB 201216C.

Detection and sky map
Fig. 1 shows the distribution of the squared angular distance,  2 , for the GRB and background events (red circles and blue squares, respectively) for the first 20 minutes of data (from  0 + 56 s to  0 + 1224 s).The significance of the VHE gamma-ray signal from GRB 201216C is 6.0 , following the prescription of Li & Ma (1983), confirming the significant detection of the GRB.For the computation of the significance, we apply cuts on  2 and hadronness.The former is the squared angular distance between the reconstructed direction of the events and the nominal position of the source, taken from Swift-BAT for GRB 201216C.The latter is a parameter which discriminates between gamma-like and background-like events, with gamma rays having hadronness values close to zero.The cuts on  2 and hadronness were optimized for a source with an intrinsic powerlaw spectrum with index  = −2, later corrected considering the absorption by the extragalactic background light (EBL) according to the model by Domínguez et al. (2011), hereafter D11.For the signal significance evaluation, the intrinsic spectral index for the cut optimization was chosen to be similar to the one found in the other GRBs detected at VHE, so without any prior knowledge of the actual value for this specific GRB (see Section 4.2).The corresponding energy threshold of the optimized cuts is 80 GeV defined by the peak of the energy distribution of the surviving simulated events.
Fig. 2 shows the test-statistics map in sky coordinates for the first 20 minutes of data.The same event cuts as for Fig. 1 are used.Our test statistic is Li & Ma (1983) equation 17, applied on a smoothed and modeled background estimation.Its null hypothesis distribution mostly resembles a Gaussian function, but in general can have a somewhat different shape or width.In the sky map, the peak position around the center is consistent with the one reported by Swift-XRT within the statistical error.The peak significance is above 6 , which corroborates the detection.

Average spectrum
The average spectrum for the first 20 minutes of observation is shown in Fig. 3.The data points are the result of an unfolding procedure following the prescription of the Bertero method described in Albert et al. (2007).The best fit to the points is instead provided by the forward folding method (Piron et al. 2001).For the event cuts optimization, the adopted spectrum is an intrinsic power-law spectrum with an index  = −3, which is close to the final estimated value (see below), later attenuated by EBL assuming the model D11 and  = 1.1.Because of the strong EBL absorption, the observed  denoted by white and blue filled points respectively.The highest energy bin is a 2 upper limit in each spectrum.The solid black and dashed grey lines represent the forward folding fits to the data points.The solid grey line is obtained from the intrinsic spectrum fit (black solid line) after the absorption by the EBL is taken into account, using the D11 model.spectrum has a steep power-law index of −5.32 ± 0.53 (stat.only) above 50 GeV.The intrinsic (EBL-corrected) spectrum is consistent with a simple power-law function and shows a harder index of −3.15 ± 0.70 (stat.only).The normalization factor at 100 GeV is (2.03 ± 0.39) × 10 −8 TeV −1 cm −2 s −1 (stat.only).The highest energy bin around 200 GeV is a 2 upper limit due to a large relative flux error about 100%.
The obtained spectrum suffers from systematic uncertainties coming from different sources.For such a steep observed spectrum, the uncertainty of the energy scale significantly affects the computed fluxes.We estimated the flux variation by shifting the light scale in the simulations during the forward folding procedure assuming the EBL model D11.We adopted a ±15% shift as prescribed in Aleksić et al. (2016).The results are shown in Table 1.When the energy scale is shifted by -15%, the observed spectrum is shifted to the low-energy side resulting in a lower flux.The spectral index of the intrinsic spectrum is softened due to the smaller attenuation by EBL at lower energies.In case of the +15% shift, the flux and the spectral index are shifted in the opposite direction.ranges from -3.19 in the -15% case to -2.17 in the +15% case.The normalization factor instead varies by a factor of 3. The spectral uncertainty originating from the energy scale is therefore significantly larger than the statistical errors.
The VHE flux of GRB 201216C is also affected by the choice between available EBL models.At such high redhift  = 1.1, EBL models show large differences in predicted attenuation factors.We compared the spectra calculated with four EBL models including D11 with the same unfolding method as the one used for Fig. 3.The three models besides D11 are Franceschini, A. et al. (2008), Finke et al. (2010), and Gilmore et al. (2012) (hereafter F08, FI10, and G12, respectively).The results are shown in Table 1.The power-law index ranges from -3.19 in the F08 case to -2.45 in the G12 case, and the normalization factor varies by a factor of 2. Also in this case, the systematic uncertainty on the parameters due to the EBL models is larger than or equal to the statistical errors.
At  = 1.1, D11 and F08 have similar attenuation values below 200 GeV, which is the maximum energy in our analysis.The attenuation discrepancy between D11 and G12 is a factor of 2 at 100 GeV and a factor of 5 at 200 GeV.Thus, the intrinsic spectrum has a larger normalization and it is harder in the G12 case than in the D11 and F08 case, as seen in Table 1.

Light Curve
The VHE energy-flux light curve between 70 GeV and 200 GeV is shown in Fig. 4. The energy flux of each time bin is obtained by integrating the EBL-corrected forward-folded spectrum with the D11 model, so that the spectral variability with time is taken into account.For each time bin, the event cut is based on the signal survival fraction of simulated events in order to increase the statistics in such short time bins.The corresponding energy threshold is around 70 GeV for all the time bins.The light curve is compatible with a power-law decay.The best fit decay index until the 5th bin excluding upper limits is −0.62 ± 0.04.
Upper limits are calculated for bins where relative flux errors are larger than 50% using the method described in (Rolke & López and from  0 +20.5 h to  0 +24.6 h (second night).Upper limits are calculated as 95% confidence level for the bins with relative errors >50%.
2001).The excess count upper limit of 95% confidence interval is calculated for each of such bins and converted into the energy flux unit by assuming the power-law spectrum with an index of -3 attenuated with the D11 model.
The systematic uncertainties considered in Sec.4.2 also affect the flux points in the light curve to a similar extent.However, since the spectral shape is not expected to change significantly during the short period of each bin of the light curve, the relative flux error is similar among all the bins.Therefore, the temporal decay index should be independent of the uncertainties as long as the spectrum is assumed to be stable.In fact, we could not detect any significant spectral changes larger than the statistical error during the time interval where the light-curve was produced.
From the analysis of the data on the second night, which spans from  0 +20.5 h to  0 +24.6 h, we found no significant excess around the position of the GRB with both the cut used in Sec.4.1 and a conventional cut optimized for the Crab Nebula.
We calculated the flux upper limit on the second night assuming an intrinsic power-law spectrum with an index of -3 and the EBL model D11.The event cut applied is the same one as used for the light curve on the first night.The upper limit of the EBL-corrected flux is shown in Fig. 4. We note that the VHE luminosity of GRB 201216C implied by MAGIC observations is fainter (a factor 10-30) than the luminosity predicted by Zhang et al. (2023) on the basis of their afterglow modeling at lower frequencies.

MULTI-WAVELENGTH DATA FROM RADIO TO GAMMA-RAY
In this Section we give an overview of the data at lower energies collected from the literature or analysed in this work, and later used for modeling and interpreting the overall emission (see Section 6).
The data is shown are Figs.5 and 6.

Radio observations
We collected radio observations from Rhodes et al. (2022).These late-time observations have been performed with e-MERLIN, the VLA, and MeerKAT, and cover the ∼ 1 − 10 GHz frequency range.There are no simultaneous detections available at higher frequencies at the time of radio detections, which span the temporal window ∼ 5 − 56 days after  0 .Rhodes et al. (2022) argue that the emission detected in the radio band is dominated by a different component as compared to the emission detected at earlier times in the optical band and in X-rays, and they suggest radiation from the cocoon as possible explanation.In our analysis we also find that radio data cannot be easily explained as synchrotron radiation from the forward shock driven by the relativistic jet, see the discussion in Section 6.We nevertheless include radio data in our analysis (star symbols in Figs. 5 and 6), verifying that the estimation of the synchrotron flux from the jet given by the modeling lies below the observed radio emission.

Optical observations: the Liverpool Telescope and VLT
The 2-m fully robotic Liverpool Telescope (LT) autonomously reacted (Guidorzi et al. 2006) to the Swift-BAT alert, and started observations from about 178 s after the burst with the IO:O5 optical camera in the SDSS-r band (Shrestha et al. 2020).The light curve of GRB 201216C optical counterpart initially displayed a flat behaviour (see Fig. 5) followed by a steepening, as revealed by VLT data gathered at 2.2 hours post-burst (Izzo et al. 2020) and by the non-detection of the afterglow in deeper LT observations at 1 day post-burst.We note that the LT photometry data was calibrated using a common set of stars present in the field of view selected from the APASS catalog.

X-ray observations
Swift-XRT started to collect data on GRB 201216C only 2966.8 s after the burst onset due to an observing constraint (Beardmore et al. 2020).Observations continued up to  0 + 4325.4 s.The unabsorbed X-ray flux integrated in the 0.3-10 keV energy range is shown in Fig. 5 (blue data points).At around 0.1 days XRT and optical data are simultaneously available and we built the spectral energy distribution (SED) around this time (Fig. 6).The XRT spectrum has been derived by analysing data between 8900 and 9300 s with the XSPEC software.Source and background spectra have been built using the automatic analysis tool6 .We modeled the spectrum with an absorbed power-law accounting both for Galactic and intrinsic metal absorption using the XSPEC models  and , respectively.The Galactic contribution is fixed to the value  H,G = 5.04 × 10 20 cm −2 (Willingale et al. 2013), while the column density in the host galaxy is a free parameter.We find that the best fit photon index is −1.67 ± 0.19 and the intrinsic column density is  H = (1.48±0.52)×10 22 cm −2 .The spectral data, rebinned for plotting purposes and de-absorbed for both Galactic and intrinsic absorption, are shown in Fig. 6 (black crosses).

Gamma-ray observations by Fermi-LAT
Fermi-LAT observations started from  0 + 3500 s and continued until the GRB position was no longer visible ( 0 + 5500 s).No signal is detected during this time window.Assuming a photon index  = −2, the estimated upper limit in the energy range 0.1-1 GeV is 3 × 10 −10 erg cm −2 s −1 (Bissaldi et al. 2020).This upper limit is included in our analysis (orange arrow in Fig. 5).
All the light curves at different frequencies are shown in Fig. 5.The Swift-BAT prompt emission light curve is also included in the figure (grey data points).The BAT flux is integrated in the 15-50 keV energy range and points are rebinned using a signal-to-noise ratio (SNR) criterion equal to seven7 .The vertical colored stripes mark the times where SEDs are built.The SEDs are shown in Fig. 6, where the MAGIC spectrum integrated between 56 s and 1224 s is also shown.

MODELING
In this Section, we discuss the origin of the emission detected by MAGIC and its connection to the afterglow emission at lower energies, from radio to X-rays.In particular, we test an SSC scenario from electrons accelerated at the forward shock.We consider a relativistic jet with initial Lorentz factor Γ 0 ≫ 1, opening angle  jet , and a top-hat geometry.The (isotropic equivalent) kinetic energy of the jet  k is related to  iso, ∼ 6 × 10 53 erg (see Sec. 2) through the efficiency for production of prompt radiation   : The details of the equations adopted to describe the dynamics, the particle acceleration, and the radiative output can be found in Miceli & Nava (2022) and are also reported in Appendix A. We summarize here the general model and the main assumptions.
The jet is expanding in an ambient medium characterised by a density described by a power-law function () ∝  − .We consider the density to be either constant ( = 0 and () =  0 ) or shaped by the progenitor stellar wind: () =  −2 ( = 2), where  is related to the mass loss rate of the progenitor's star  and to the velocity of the wind  w by  = /4   p   ( p is the mass of the proton).We normalize the value of  to a mass loss rate of 10 −5 solar masses per year and a wind velocity of 10 3 km s −1 :  = 3 × 10 35  ★ cm −1 .We assume that ambient electrons are accelerated at the forward shock into a power-law distribution / ∝  −  from  min and  max .The bulk Lorentz factor of the fluid just behind the shock is assumed to be constant (Γ = Γ 0 ) before the deceleration and described by the solution given by Blandford & McKee (1976) (Γ = Γ BM ) during the deceleration (note that the equation given by Blandford & McKee (1976) describes the Lorentz factor of the shock Γ sh , which we relate to the Lorentz factor of the fluid using Γ = Γ sh / √ 2).The two regimes are smoothly connected to obtain the description of the bulk Lorentz factor of the fluid just behind the shock as a function of shock radius.
To infer the particle distribution and the photon spectrum at any time  we numerically evolve the equations describing the electron and photon populations including adiabatic losses, synchrotron emission and self-absorption, Inverse Compton emission and  −  annihilation and pair production.To relate the comoving properties computed by the code to the observed one, we assume that the emission received at a given observer time is dominated by electrons moving at an angle cos  =  from the line of sight to the observer, where  is the velocity of the shocked fluid.
Before presenting the results of the numerical modeling, we discuss some general considerations that can be inferred using analytic approximations from Granot & Sari (2002).Fig. 5 shows that the optical flux is nearly constant up to at least 5 × 10 −3 d.At later times this behaviour breaks into a steeper temporal decay.We take as reference value for the break time ∼ 10 −2 d.This behavior of the optical light curve can be explained if the break frequency  m (i.e. the typical photon energy emitted by electrons with Lorentz factor  min ) is crossing the  band.The nearly constant flux before the crossing time is indicative of a wind-shaped external medium (i.e.,  = 2).The preference for a wind-like medium is also supported by the lack of a phase of increasing flux in the MAGIC observations, which start as early as ∼60 s after the onset of the prompt emission.Since  m ∝  −1.5 , we expect  m ∼ 1 GHz at ∼ 50 d.Observations at 1.3 GHz do not allow to constrain the peak time, but we notice that they are consistent with the presence of a peak around 50 days (a zoom on radio observations can be found in figure 1 of Rhodes et al. 2022 and it shows that, considering the errors, the flux is consistent with being constant at about 50 days).The radio SED at this time (see Fig. 6 and also Rhodes et al. 2022) shows that the self-absorption frequency  sa must be below 1 GHz.Since  sa ∝  −3/5 , this implies that during the time spanned by observations  sa <  m .The fast increase of the 1.3 GHz flux (with temporal index ≳ 5, as reported in Rhodes et al. 2022) however implies that observations at this time are below the self-absorption frequency (otherwise the flux at 1.3 GHz should be constant), constraining  sa (50 d) ∼1 GHz.
To summarise, the scenario implied by optical and radio observations invokes a jet expanding in a wind-like density and producing a synchrotron spectrum with  sa <  m , and  m crossing the optical  band at ∼ 10 −2 d and the 1 GHz frequency at ∼ 50 d.We now check the consistency of this interpretation with X-ray observations.Imposing  m (54 d) =  sa (54 d) = 1 GHz and the flux  ( sa , 54 d) = 2 × 10 −18 erg cm −2 s −1 , and using equations for the break frequencies and flux in a wind-like medium from Granot & Sari (2002), it is possible to derive the values of  k ,  B , and  ★ as a function of  e , for fixed values of .For  = 2.2 we find  k,52 ≃  e ,  B ≃ 2.3 × 10 −6  −5 e and  ★ ≃ 8.8  2 e .This shows that the requirement that the   spectrum peaks at  sa = 1 GHz at 54 days limits the energy to a low value  k < 10 52 erg, inconsistent with the large flux detected in the X-ray afterglow.In particular we find that the X-ray band is always above the cooling frequency  c for different assumptions on  e .Moreover, in this range, the predicted flux is at least one order of magnitude below the detected X-ray flux.This statement is quite robust, as the flux in this band weakly depends on  B , does not depend on  ★ , and is proportional to  k  e .Pushing  e to large values (close to one) improves the situation, at the expense of a very small  B , implying a large SSC component.This solution is ruled out by MAGIC observations.
Being unable to find a scenario that explains all the available data as synchrotron and SSC emission from the forward shock driven by a relativistic jet, we consider the possibility that late-time radio emission is dominated by a different component, as also concluded by Rhodes et al. (2022), which identify in a wider mildly (or non) relativistic cocoon the origin of the radio emission.We then restrict the modeling to the MAGIC, X-ray and optical data, requiring that the flux at 1-10 GHz from the narrow relativistic jet is below the observed flux.
We performed numerical calculations of the expected synchrotron and SSC radiation and their evolution in time for wide ranges of values of the parameters  k ,  e ,  B , , (),  jet , and Γ 0 .The investigated range of values for each parameter is reported in Table 2.The numerical calculations confirm the considerations derived from analytic estimates.In particular, we neither find a solution for a homogeneous medium nor for a complete description of radio to GeV observations.Assuming a wind-like density profile, we find that the observations can be well described as synchrotron and SSC radiation.In particular, once the request to model also radio observations with forward shock emission from the relativistic jet is abandoned, the X-ray flux can be explained by increasing the assumed value of the jet energy, which also moves the self-absorption frequency to lower energies.An example of modeling is provided in Figs. 5 and 6, where observations (corrected for absorption in the optical and X-ray band) are compared to the light curves and spectra predicted with the following parameters:  k = 4 × 10 53 erg,  e = 0.08,  B = 2.5 × 10 −3 ,  ★ = 2.5 × 10 −2 ,  = 2.1, Γ 0 = 180, and  jet = 1 • the values are listed also in Table 2.The jet opening angle is broadly constrained by the need to not overproduce the radio flux.The inferred value points to a narrow jet, with opening angle in the low-value tail of distributions of inferred jet opening angles for long GRBs (Chen et al. 2020).We note that a similarly small ( jet ∼ 0.8 • ) value for the jet opening angle has been inferred for the TeV GRB 221009A LHAASO Collaboration et al. (2023).The inferred jet kinetic energy implies an efficiency of the prompt emission   ≃ 60%.In agreement with the steep optical spectrum reported by Vielfaure et al. (2020), this model implies an extinction of 4.6 magnitudes in the  ′ band, which is well in excess of the Galactic contribution ( ( −) = 0.05).As it can be seen from Fig. 5, the onset of the deceleration occurs at  obs ≲ 200 s, where the X-ray and TeV theoretical light curves steepen from an almost flat to a decaying flux.
The steepening of the optical light curve instead occurs at ∼ 10 3 s because, as already commented, it is determined by the  m frequency crossing the  band.In this interpretation, the frequency  m is initially above the optical band (see the brown SED in Fig. 6) and then moves to lower frequencies crossing the optical and explaining the steepening in the light curve.X-ray observations lie just above the cooling frequency, but the X-ray spectrum remains harder than expected due to the role of the Klein-Nishina cross section.We also computed the expected SED averaged between 56 s and 1224 s, where the MAGIC spectrum (see Fig. 3) is computed.The model SED is reported in Fig. 6 (green curve, to be compared with the MAGIC data, green circles).We find that the  −  internal absorption plays a minor role in shaping the spectrum: the flux reduction at 200 GeV is about 25%.In the same figure it is also possible to see the expected location of the maximum energy of synchrotron photons, initially located at 10 GeV at the time of the first SED, and then moving towards lower energies.Assuming diffusive shock acceleration proceeding at the maximum rate rules out a synchrotron origin for the photons detected by MAGIC.Rhodes et al. 2022).The XRT spectral data points estimated around 9000 s are also shown.Green circles show the MAGIC spectrum averaged between 56 and 1224 s (Fig. 3).The theoretical SED to be compared with the MAGIC spectrum is the green curve, which shows the predicted spectrum (synchrotron + SSC) averaged in the same time window (56-1224 s).

CONCLUSIONS
In this paper, MAGIC analysis results on GRB 201216C and their interpretation were presented.The GRB afterglow was observed at early times (∼ 10 2 − 10 3 s) by MAGIC for a total of ∼2.5 h during the first night and detected at the level of 6 in the first 20 minutes.This is the second firm detection of a GRB with the MAGIC telescopes after GRB 190114C, and also the farthest VHE source detected to date.Both the observed and intrinsic average spectra can be well described by a power-law.A time-resolved analysis was also performed, in order to evaluate the temporal behavior of VHE emission.The obtained light curve shows a monotonic power-law decay, indicating a probable afterglow origin of the VHE emission.Multi-wavelength data were also collected by other ground and space-based instruments.Unfortunately, most of them are not contemporaneous to the first MAGIC observation time window.In other cases, as for Fermi-LAT, the GRB could not be detected.Like other GRBs detected in the VHE range, multi-wavelength data was used to perform a modeling of the broadband emission.In this manuscript a synchrotron and SSC radiation model at the forward shock in the afterglow was considered.SEDs built at different times show that synchrotron photons can reach a maximum energy of 10 GeV about three minutes after the GRB onset.The emission detected by MAGIC reaches higher energies, and can therefore be explained by the SSC component of the model.By comparing analytic estimates and the numerical modeling, evidence for the need of a different component at the origin of late-time radio emission is found, in agreement with previously published studies on this GRB.Both observations and modeling support a wind-like medium, as expected in the case of a long GRB.The best fit model parameters are found to be consistent with those estimated in previous studies of GRB afterglows without VHE detection.This proves the flexibility of the SSC scenario in describing the VHE emission of GRBs.
Like other VHE detected GRBs (GRB 180720B and GRB 190114C), 201216C was a bright GRB, allowing for a detection in spite of the high redshift.Once again, the rapid response and low-energy threshold of the MAGIC telescopes to GRB alerts was crucial to detect the VHE emission in the early afterglow phase.Altogether, the detection by MAGIC and other experiments of several bursts so far suggests that VHE emission is common both in high and low-luminosity GRBs.Other VHE detected GRBs showed a correlation between the intrinsic emission in the X-ray and VHE bands, where a similar time decay and flux value were observed.In the case of GRB 201216C such a direct comparison cannot be performed given the lack of contemporaneous data in the two bands.The extrapolation of the X-ray flux into the first MAGIC time window, assuming a smooth power-law behavior typical of the afterglow phase, shows that the VHE flux is lower than the X-ray one.However, one should take into account the rather narrow energy range of the VHE detection due to the large absorption caused by the EBL.Number JP21K20368.L. Nava acknowledges partial support from the INAF Mini-grant 'Shock acceleration in Gamma Ray Bursts'.The Liverpool Telescope is operated on the island of La Palma by Liverpool John Moores University in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias with financial support from the UK Science and Technology Facilities Council (STFC) under UKRI grant ST/T00147X/1.Manisha Shrestha and Iain Steele thank UKRI/STFC for financial support (ST/R000484/1).We would like to thank the STARGATE collaboration for the confirmation of the redshift of the source via private communication.

Figure 1 .
Figure 1. 2 distribution for the first 20 minutes of observation, see the main text for the definition.Both GRB (red circles) and background events (blue squares) are shown.The vertical dashed black line shows the value of the cut in  2 used for the calculation of the significance.

Figure 2 .
Figure 2. Test-statistics sky map for the first 20 minutes of observation.The cross marker shows the position of GRB 201216C reported by Swift-XRT.The white circle shows the MAGIC point spread function corresponding to 68% containment.

Figure 3 .
Figure3.Observed and EBL-corrected spectra for GRB 201216C as measured by the MAGIC telescopes during the first 20 minutes of observations, denoted by white and blue filled points respectively.The highest energy bin is a 2 upper limit in each spectrum.The solid black and dashed grey lines represent the forward folding fits to the data points.The solid grey line is obtained from the intrinsic spectrum fit (black solid line) after the absorption by the EBL is taken into account, using the D11 model.

Figure 4 .
Figure 4. EBL-corrected energy-flux light curve between 70 GeV and 200 GeV from  0 +56 s to  0 +40 min (first night, divided into five time bins)and from  0 +20.5 h to  0 +24.6 h (second night).Upper limits are calculated as 95% confidence level for the bins with relative errors >50%.

Figure 5 .
Figure 5. Multi-wavelength light curves of GRB 201216C.Both the X-ray and optical observations have been corrected accounting for absorption.MAGIC data points are EBL-corrected.Upside down triangles represent upper limits.Solid curves show the best fit model obtained in a synchrotron -SSC forward shock scenario.Different colors refer to the different wavelengths where observations are available (see the legend).The modeling is obtained with the following parameters:  k = 4 × 10 53 erg,  e = 0.08,  B = 2.5 × 10 −3 ,  ★ = 2.5 × 10 −2 ,  = 2.1, Γ 0 = 180, and  jet = 1 • .Vertical lines mark the times where SED have been built (see Fig. 6).

Figure 6 .
Figure6.SEDs of GRB 201216C at different times.Different colors for curves and data points refer to different times (see the legend).The times where the SEDs are calculated are also marked in Fig.5with vertical stripes.Solid curves show the synchrotron and SSC theoretical spectra for the same parameters used for Fig.5.De-absorbed optical data in the  ′ filter are marked with square symbols, while star symbols are observations at 1.3 GHz and 10 GHz at 54.5 d and 53 d, respectively (fromRhodes et al. 2022).The XRT spectral data points estimated around 9000 s are also shown.Green circles show the MAGIC spectrum averaged between 56 and 1224 s (Fig.3).The theoretical SED to be compared with the MAGIC spectrum is the green curve, which shows the predicted spectrum (synchrotron + SSC) averaged in the same time window (56-1224 s).

Table 1 .
The obtained power-law index Fitted power-law spectral parameters of the 20-minute average spectrum using different scales of the Cherenkov light amount and different EBL models.The tested EBL models are D11, F08, FI10, and G12 with the nominal light scale.The tested light scales are nominal, -15%, +15% with the D11 EBL model.The normalization energy is fixed to 100 GeV.The errors are statistical only.The resulting systematic errors are reported in the main text.

Table 2 .
List of the input parameters for the afterglow model.For each parameter, the range of values investigated by means of the numerical model are listed in the second column.Solutions are not found for an homogeneous density medium ( = 0).The last column list the values that better fit the observations and used to produce the model light-curves and model SEDs in Figs. 5 and 6.