Prompt emission and early optical afterglow of VHE detected GRB 201015A and GRB 201216C: onset of the external forward shock

For GRB 201015A, Swift could not slew until T₀ +51.6 minutes due to an observing constraint (D’Elia et al., 2020). The Swift X-ray telescope (XRT, Burrows et al. 2005) began observations of GRB 201015A at 23:43:47.2 UT on 2020 October 15, $\sim$ 3214.1 s post burst. Swift-XRT detected an uncatalogued fading X-ray source at the following location: RA, Dec = 354.32067, +53.41460 degrees (J2000) with an uncertainty radius of 3.8^" (Kennea et al., 2020). Swift-XRT observed this source up to $\sim{1.8}$ $\times~{}10^{3}$ ks after the BAT detection. All observations were obtained in Photon Counting (PC) mode.

For GRB 201216C, Swift-XRT began observations at 23:56:58.5 UT on 2020 Dec 16, $\sim$ 2966.8 s post burst. Swift-XRT detected a new fading X-ray source at the following location: RA, Dec = 16.37114, +16.51659 degrees (J2000) with an uncertainty circle of 3.5^" (Campana et al., 2020). Swift-XRT observed this source up to 2.2 $\times~{}10^{3}$ ks after the initial detection. Window timing (WT) mode was used for the first $\sim$ 25 s of observations, and the remaining observations were obtained in PC mode.

X-ray afterglow data analysis: In this work, we obtained the Swift-XRT data products from the XRT repository provided by the University of Leicester (Evans et al., 2007, 2009). We modeled the X-ray afterglow light curves of both the bursts using power-law and broken power-law empirical models to constrain their decay rates. On the other hand, to constrain the spectral indices, we modeled the XRT spectra of both the bursts in the energy range of 0.3-10  keV using the X-ray Spectral Fitting Package (XSPEC; Arnaud, 1996). We fixed the Galactic hydrogen column density to be $\rm NH_{\rm Gal}$= $3.60\times 10^{21}$, and $5.04\times 10^{20}$ cm^-2 for GRB 201015A and GRB 201216C, respectively (Willingale et al., 2013). A more detailed method for the standard temporal and spectral XRT analysis is given in Gupta et al. (2021b).

GRB 201015A: The optical afterglow of GRB 201015A was first reported by the MASTER robotic telescope at the position RA= 23:37:16.42, DEC= +53:24:55.8 (Lipunov et al., 2020). In this paper, we present our optical observations from the 25cm FRAM-ORM (Jelinek et al., 2020a), Burst Observer and Optical Transient Exploring System (BOOTES) (Hu et al., 2020), and 3.6m Devasthal Optical Telescope (DOT), along with additional optical data from the GCN circulars. Our photometric observations are listed in Table LABEL:table:15A_optical of the appendix.

BOOTES-network: The BOOTES-network followed GRB 201015A with three robotic telescopes at BOOTES-1 (INTA-CEDEA) station in Mazagon (Huelva, Spain) and BOOTES-2/TELMA station in La Mayora (Malaga, Spain). The BOOTES-1B performed one epoch of early observations at 22:50:48 on 2020 October 15, and the afterglow is clearly seen in the images. The BOOTES-1A performed two epoch observations: the first was the quick follow-up to the trigger and the second was the late follow-up on the next day. The BOOTES-2/TELMA also followed this event on 2020-10-16T22:15:13. The afterglow was not visible in both BOOTES-1A’s and BOOTES-2’s images.

A series of images were obtained with BOOTES-1B robotic telescope in the clear filter with exposures of 1 s, 10 s, and 60 s for GRB 201015A. The image preprocessing (bias and dark-subtracted, flat-fielded, and cosmic-ray removal) was done using custom IRAF routines. The photometry was carried out using the standard IRAF package. The images were calibrated with nearby comparison stars from the Pan-STARRS catalog, which were imputed to R-band through the transformation equation from Lupton (2005)⁹⁹9http://www.sdss.org/dr12/algoritms/sdssUBVRITransform/ for the data in the clear filter. The obtained magnitudes are listed in Table LABEL:table:15A_optical of the appendix.

FRAM-ORM: The FRAM-ORM is a 25cm f/6.3 robotic telescope facilitated with B, V, R, and z filters and a custom Moravian Instruments G2-1000B1 camera based on a back-illuminated CCD47-10 chip. We carried out the earliest optical observations of GRB 201015A using the 25cm FRAM-ORM robotic telescope located at La Palma, Spain. We obtained a series of frames with exposure times of 20 s each in the clear filter, beginning on 2020 October 15 at 22:50:50.8 UT (37.6 s after the BAT trigger). We have clearly identified the source mentioned by Lipunov et al. (2020). A finding chart of FRAM-ORM observation of GRB 201015A is given in Figure 1.
3.6m DOT: We performed the optical observations of the afterglow of GRB 201015A starting on 2020 October 16 at 13:09:07.2 UT ($\sim$ 0.6 days post burst) using 3.6m DOT located at the Devasthal observatory, which is part of the Aryabhatta Research Institute of Observational Sciences (ARIES), India. We acquired multiple frames in B, V, R, and I filters (with an exposure time of 10 minutes in each) using the 4K $\times$ 4K CCD IMAGER (Kumar et al., 2022b) mounted on the main port of 3.6m DOT. In stacked DOT images, we clearly detect the optical afterglow of GRB 201015A consistent with the error region of NOT observations Malesani et al. (2020). A finding chart of 3.6m DOT observation of GRB 201015A is given in Figure 15 of the appendix. We performed the DOT image reduction using the IRAF package. We first applied zero correction and flat fielding to the raw images taken from the telescope. After the removal of cosmic-ray hits, we stack the images to create a single image. We use the IRAF package to perform the aperture photometry. For the photometric calibration of GRB 201015A, standard stars in the Landolt standard field PG 0231 were observed along with the GRB field in UBVRI bands. The R-band finding chart of GRB 201015A and the secondary stars marked with S1-S14 are shown in the Figure 15. The calibrated magnitudes of secondary stars are listed in Table 10 of the appendix, and the calibrated magnitudes of GRB 201015A are listed in Table LABEL:table:15A_optical of the appendix.

GRB 201216C: In the case of GRB 201216C, Izzo et al. (2020) carried out the optical follow-up observations in the Sloan g${}^{{}^{\prime}}$, r${}^{{}^{\prime}}$, z${}^{{}^{\prime}}$ filters using VLT telescope, starting at 01:18:47 UT on 2020 December 17. They first reported the detection of an uncatalogued optical source at location RA: 01:05:28.980, Dec: +16:31:00.0 (J2000.0), consistent with the Swift-XRT enhanced location (Campana et al., 2020).

FRAM-ORM: We performed the earliest optical follow-up observations to the alert of GRB 201216C using the FRAM-ORM telescope. We obtained a series of frames with an exposure times of 20 s each in the clear filter, starting on 2020 December 16 at 23:08:04.3 UT (31.6 s after the BAT trigger) (Jelinek et al., 2020b). We immediately detected an optical transient consistent with the location reported by Izzo et al. (2020). A finding chart of FRAM-ORM observation of GRB 201216C is given in Figure 1. A log of photometric observations of the afterglow of GRB 201216C is presented in Table 12 of the appendix.

The observed frames for both the bursts have been processed through a difference imaging pipeline based on the HOTPANTS image subtraction code to remove the influence of nearby stars on the photometric measurements and then corrected for bias, flats, and cosmic-rays. The photometry has been carried out using the DAOPHOT package. The measured magnitudes have been calibrated using field stars in the PanSTARS DR1 catalog. Our optical observations reveal an early rise in the light curves of both bursts, reaching their maximum and followed by normal decay until the end of our observations. A log of our photometric observations is given in Table LABEL:table:15A_optical of the appendix.

The Swift-BAT energy-resolved mask-weighted light curve of GRB 201015A along with HR evolution in 25-50  keV and 50-100  keV energy range is shown in the top panel of Figure 2. The Swift-BAT prompt emission light curve of GRB 201015A has a short-soft emission starting from T₀ and ending at T₀ +1 s, followed by a weak-soft tail emission that lasts till $\sim$ T₀ +10 s (Markwardt et al., 2020). The time-integrated spectrum (T₀-0.21 to T₀ +11.57 s) of GRB 201015A is explained using the simple power-law model with power index = -2.43 $\pm$ 0.25 (see Figure 8 of the appendix).

The background-subtracted energy-resolved Fermi-GBM light curve of GRB 201216C along with the hardness-ratio (HR) evolution is shown in the bottom panel of Figure 2. The Fermi-GBM prompt emission light curve consists of a broad, structured peak with a T₉₀ duration of $\sim$ 29.9 s (in 50-300  keV). Similar to the Fermi-GBM, Swift-BAT observations also reveal multiple peaks in the mask-weighted prompt light curve. According to BAT observations, the energy fluence of the burst is 4.5 $\pm 0.1\times 10^{-5}$ erg $\rm cm^{-2}$ in the 15-150  keV energy range (Ukwatta et al., 2020). Figure 2 also shows the time interval used for the time-averaged spectral analysis. The time-averaged Fermi-GBM spectrum (from T₀-0.503 s to T₀ +47.09 s) is best modeled using a Band+BB function (see Figure 8 of the appendix) and the best fit spectral parameters are presented in Table 4 of the appendix. The burst was significantly bright: the observed fluence is (2.00 $\pm$ 0.10) $\times 10^{-4}$ erg $\rm cm^{-2}$ in the 10-10⁴  keV energy band over the time-averaged interval. This fluence value is among the top 2 percent of the brightest GRBs observed by Fermi-GBM. A comparison between the time-averaged energy spectrum of the empirical (Band, CPL) and physical (Synchrotron) models are shown in Figure 3.

GRBs are classified into two families based on T₉₀ duration, however, there is one more fundamental difference between both the classes, i.e., spectral hardness. Long GRBs are expected to be softer (lower $E_{\rm p}$ value), on the other hand, short bursts are usually harder (higher $E_{\rm p}$ value). To know the true class of GRB 201015A and GRB 201216C, we placed both bursts in the time-integrated T₉₀-spectral hardness plane along with the other long and short GRBs obtained from the Fermi-GBM and Swift-BAT catalogs.

We calculated the hardness ratio (HR) for GRB 201015A using the ratio of fluence in hard (50-100  keV) and soft (25-50  keV) energy ranges and found its value equal to 0.72. Comparing the HR of GRB 201015A with the Swift-BAT catalog, we note that GRB 201015A is one of the softest bursts ever observed by the Swift mission (see Figure 9 of the appendix). In the case of GRB 201216C, we calculated the HR using the ratio of counts in hard (50-300  keV) and soft (8-30  keV) energy ranges to be 1.11, this feature is similar to long/soft GRBs. The distribution of HR as a function of T₉₀ for both the bursts are given in Figure 9 of the appendix.

Furthermore, we also placed both the bursts in the time-integrated $E_{\rm p}$-T₉₀ distribution of long and short GRBs. As for GRB 201015A there was no onboard observation by Fermi and no public Fermi data is available for ground targeted search, we have used Swift-BAT observations of GRB 201015A to constrain the time-integrated peak energy. However, the BAT time-integrated spectra of GRBs are usually modeled by simple power-law functions due to the instrument limited spectral coverage (15-150  keV). In the case of GRB 201015A, the time-integrated spectrum is also fitted using a power-law function (see § 3.1.1). We estimated the $E_{\rm p}$ value using the known correlation between the Swift-BAT fluence and time-integrated $E_{\rm p}$ (Zhang et al., 2020), i.e., $E_{\rm p}$ = [fluence/($\rm 10^{-5}erg~{}cm^{-2}$)]^0.28 $\rm\times 117.5^{+44.7}_{-32.4}$  keV $\approx$ 41.37${}^{+15.74}_{-11.41}$  keV. The estimated time-integrated $E_{\rm p}$ value of GRB 201015A is consistent with those observed for long GRBs. In the case of GRB 201216C, we fitted the time-integrated spectrum and calculated the peak energy ($E_{\rm p}$ = 352.31${}^{+12.77}_{-12.74}$  keV) using the best fit model. The T₉₀-$E_{\rm p}$ distribution for both the GRBs is shown in Figure 9 of the appendix along with other data points obtained from Fermi-GBM catalog (Goldstein et al., 2017). We fitted the complete distribution obtained from Fermi-GBM catalog using the Bayesian Gaussian mixture model (BGMM) algorithm and calculated the probability of GRB 201216C being a long burst as 98.8 %.

The cosmological corrected time-integrated peak energy $E_{p,z}$ = (1+z)$E_{\rm p}$ of the prompt emission is correlated with the isotropic equivalent energy $E_{\rm\gamma,iso}$, and isotropic peak luminosity $L_{\rm\gamma,iso}$. The former is known as Amati correlation (Amati, 2006) and the later is known as Yonetoku correlation (Yonetoku et al., 2004). In the case of GRB 201015A, we have used equation 6 of Fong et al. (2015) to calculate $E_{\rm\gamma,iso}$ due to the limitation of Swift-BAT spectral coverage. On the other hand, we have used the best fit Fermi-GBM time-integrated spectral model to calculate the $E_{\rm\gamma,iso}$ for GRB 201216C. The Amati correlation for both the GRBs along with a sample of long and short GRBs taken from Minaev & Pozanenko (2020) is presented in Figure 9 of the appendix. Similarly, we calculated the $L_{\rm\gamma,iso}$ values and placed GRB 201015A and GRB 201216C in the $E_{p,z}$-$L_{\rm\gamma,iso}$ plane, as shown in Figure 9 of the appendix. We noticed that both the bursts satisfied the Amati and Yonetoku correlations. The calculated isotropic equivalent luminosity values for both the bursts suggest that GRB 201015A is an intermediate luminous GRB; on the other hand, GRB 201216C is a luminous GRB.

Distribution of spectral parameters: The mean values and standard deviation for each spectral parameter obtained using Cutoff power-law, Band, Cutoff power-law + BB, Band+ BB and Synchrotron models are listed in Table 8 of the appendix. The mean values of the low energy spectral indices $\alpha_{\rm CPL}$ obtained using Cutoff power-law and $\alpha_{\rm pt}$ obtained using Band spectral modeling of GRB 201216C are $-1.10\pm 0.14$ and $-1.06\pm 0.13$, respectively. These values are consistent with the typical average value of the low energy spectral index ($\alpha_{\rm pt}$ $\sim~{}-1$) of GRBs (Preece et al., 2000). Similarly, the averaged value of spectral peak energy $E_{\rm p}$ and high energy photon index $\beta_{\rm pt}$ obtained using Band are $339.43\pm 119.39$ keV, and $-2.75\pm 0.33$, respectively. These values are also consistent with the typical average value of $E_{\rm p}$ and $\beta_{\rm pt}$ of GRBs. The averaged values of physical synchrotron spectral parameters for GRB 201216C are following: magnetic field (B) = $96.00\pm 45.82$ G, index (p) = $4.18\pm 0.54$, and Lorentz factor corresponds to the electron’s cooling time $\gamma_{cool}$ = $(6.08\pm 4.81)\times 10^{5}$.

Evidence for thermal component: In our time-resolved spectral analysis, we first fitted each binned spectra using empirical Band and Cutoff power-law models individually. To search for the presence of an additional thermal component in the spectrum, we added the Blackbody model along with Band and Cutoff power-law models, and calculated the difference of the DIC values. Negative values of $\rm\Delta DIC_{\rm Band}$ and $\rm\Delta DIC_{\rm CPL}$, implies that the additional Blackbody component improves the fit statistic. Furthermore, in the case of any particular bin, if $\rm\Delta DIC<-10$, it suggests a significant amount of thermal component present in that particular bin. Figure 10 of the appendix shows the evolution of $\rm\Delta DIC$ within the burst interval for Band + Blackbody model. For all the bins (except the last bin), $\rm\Delta DIC$ values are negative, indicating that the addition of thermal components improves the spectral fitting. There are 10 bins ($\sim$ 37 %) for which $\rm\Delta DIC<-10$, suggesting for presence of thermal component in the spectrum. In light of this, we suggest that GRB 201216C is a hybrid (non-thermal+ quasi-thermal) burst. The thermal components are dominating during the initial and bright phases of the burst.

The evolution of spectral parameters: The studies of spectral evolution are a very powerful tool to probe the emission process responsible for the prompt emission. The peak energy of the Band function shows four different possible types of spectral evolution. 1. Hard to soft pattern (Norris et al., 1986), 2. Soft to hard pattern (Kargatis et al., 1994), 3. flux tracking pattern (Golenetskii et al., 1983) and 4. Chaotic pattern. On the other hand, the low energy spectral index of the Band function also changes with time, however, it does not show any particular trend. There are some recent studies suggesting a flux tracking pattern of $\alpha_{\rm pt}$, supporting the double tracking trend (Li et al., 2019). However, most of the spectral evolution studies are performed for single pulsed bursts. The prompt light curve of GRB 201216C shows a more complex multi-pulsed structure, and the evolution of empirical and physical spectral parameters are very interesting. Figure 4 shows the evolution of the spectral parameters ($E_{\rm p}$, and $\alpha_{\rm pt}$) of GRB 201216C obtained using empirical Band function. The evolution of $E_{\rm p}$ of GRB 201216C shows a flux tracking trend, i.e., $E_{\rm p}$ increases and decrease as flux increases and decrease in respective bins. We noticed that the $\alpha_{\rm pt}$ values are changing with time and do not exceed the expected values of spectral indices of synchrotron fast and slow cooling cases. We have also shown the evolution of $\alpha_{\rm pt}$ and $E_{\rm p}$ together in Figure 4, and we can see that the evolution patterns are quite similar throughout the emission. Next, we study the spectral evolution of parameters obtained from the physical synchrotron modeling. Figure 5 shows spectral evolution of the magnetic field strength (B), and spectral index of electrons ($p$). We noticed that the magnetic field strength is also following the intensity of the burst. We could not confirm the evolution trend of $p$ due to large associated errors. In light of the above, we suggest that multi-pulsed GRB 201216C have flux tracking characteristics. We further investigated the correlation among these spectral parameters.

Spectral parameters correlations: Studying the correlations between the different spectral parameters obtained using the time-resolved spectral analysis using empirical and physical modeling gives an important clue about the physics of GRBs and jet composition. We calculated the correlation between Band and Band+Blackbody spectral parameters: log (Flux) and log ($E_{\rm p}$), log (Flux) and $\alpha_{\rm pt}$, log ($E_{\rm p}$) and $\alpha_{\rm pt}$, and log (Flux) and k$T$. We found a high degree (the correlation coefficient ranges from 0.60-0.79) of correlation between the Band spectral parameters, i.e., log (Flux) and log ($E_{\rm p}$), log (Flux) and $\alpha_{\rm pt}$, and log ($E_{\rm p}$) and $\alpha_{\rm pt}$. We also found a high degree of correlation between log (Flux) and k$T$ for the Band+ Blackbody model. The correlation results for the Band and Band+ Blackbody models (Pearson correlation) are listed in Table 9 and Figure 11 of the appendix. We also studied the correlation between physical spectral parameters (log (Flux)-log (B), log (Flux)-$p$, and log (B)-$p$) obtained using Synchrotron model. We found a very high degree (the correlation coefficient is $\geq$ 0.80) of correlation between log (Flux)-log (B) and medium correlation (the correlation coefficient ranges from 0.40-0.59) between log (B)-$p$, however, we found a low degree (the correlation coefficient ranges from 0.20-0.39) of correlation between log (Flux)-$p$. The correlation results obtained using physical modeling are shown in Figure 12 and in Table 9 of the appendix.

In addition, we also studied the correlation between empirical and physical models parameters: log (B)-log ($E_{\rm p}$), log (B)-$\alpha_{\rm pt}$, and $p$-$\alpha_{\rm pt}$. The correlation results are shown in Figure 13 and Table 9 of the appendix. We found a very high degree of correlation between log (B) of physical synchrotron model and log ($E_{\rm p}$) of empirical Band function. We also found a medium degree of correlation between log (B) and $\alpha_{\rm pt}$, however, there is a low degree of correlation between $p$ and $\alpha_{\rm pt}$.

The X-ray afterglow light curves of GRB 201015A and GRB 201216C do not show any flare, bump, break, or plateau-like activities (see Figure 6).
In the case of GRB 201015A, we initially tried to fit the X-ray light curve using a simple power-law model (temporal index = -2.36${}^{+0.17}_{-0.26}$). The calculated XRT temporal index is consistent with the temporal index reported on the XRT catalogue page using count rate light curve fitting¹⁰¹⁰10https://www.swift.ac.uk/xrt_live_cat/01000452/. However, the model shows a significant deviation from the last observed data point (chi/dof =1.34). We noted that the last two data points are considered unreliable because no centroid could be determined. We calculated the XRT temporal decay index = -2.27${}^{+0.43}_{-0.58}$ excluding these points. We found that both the calculated XRT decay indices are consistent, considering the large associated errors. We also used the broken power-law function to fit the light curve and obtained chi/dof = 0.27, indicating that the model is over fitting the data. Further, we used F-test to find the best fit model among the two empirical models. The F-test suggests that the simple power-law model is better fitting the X-ray afterglow light curve of GRB 201015A, consistent with those reported on the XRT catalogue page. The X-ray light curve and the best fit model are shown in Figure 6. For the XRT spectral analysis, we fitted the late time spectrum (T₀ +4000 to T₀ +22019 s) to constrain the intrinsic hydrogen column density ($\rm NH_{\rm z}$) of GRB 201015A and calculated $\rm NH_{\rm z}$ = $5.56\times 10^{21}$ cm^-2. The time-averaged X-ray spectrum (T₀ +3300 to T₀ +1800000 s) of GRB 201015A is described using an absorbed power-law model with a spectral index = 1.10${}^{+0.25}_{-0.25}$. Further, we divided the XRT light curve into three segments based on available observations. For the first (T₀ +3300 to T₀ +4800 s), second (T₀ +10000 to T₀ +15000 s), third (T₀ +10000 to T₀ +1800000 s) time bins, we obtained the spectral indices = 1.07${}^{+0.31}_{-0.30}$, 1.10${}^{+0.53}_{-0.55}$, and 1.38${}^{+0.49}_{-0.47}$, respectively. Our analysis does not find any noticeable evolution among the spectral indices within the observed duration considering large values of associated errors.
The X-ray light curve of GRB 201216C is shown in Figure 6. The light curve (@ 10  keV) decay as a power-law with a decay index of $\alpha_{x}$ = $-2.21_{-0.11}^{+0.10}$ with (chi/dof =3.35). For the XRT spectral analysis, we calculated the $\rm NH_{\rm z}$ of the host using the late time spectral fitting (T₀ +4001 to T₀ +37346 s) and found $\rm NH_{\rm z}$ = $2.71\times 10^{22}$ cm^-2. The joint PC and WT mode time-averaged X-ray spectrum (T₀ +2900 to T₀ +17000 s) of GRB 201216C is described using an absorbed power-law model with a spectral index = 0.97${}^{+0.05}_{-0.05}$. Further, we created the spectrum for individually segmented WT and PC mode observations using the Swift Science Data Centre web-page¹¹¹¹11https://www.swift.ac.uk/xrt_spectra/addspec.php?targ=01013243&origin=GRB and performed the spectral fitting. For the WT (T₀ +2900 to T₀ +3000 s), and PC (T₀ +3300 to T₀ +17000 s) time bins, we obtained the spectral indices = 0.87${}^{+0.22}_{-0.22}$, and 0.90${}^{+0.10}_{-0.09}$, respectively. Our analysis indicates that the spectral index is not changing with time (see Figure 6).

For both these VHE detected bursts, we can examine the early afterglow behavior using the earliest optical observations from the 0.25m robotic FRAM-ORM and BOOTES telescopes (Jelinek et al., 2020a, b; Hu et al., 2020).

The optical light curve of GRB 201015A is highly rich in features. Following an early decay, the light curve has a smooth optical bump, which may be either due to reverse shock emission or the onset of afterglow in the forward shock. We calculated the decay rate for the very early optical emission $\alpha_{\rm o1}=-0.68\pm 0.15$. To characterize the nature of the early bump, we fitted the optical bump using a smoothly broken power-law function, given in equation 1.

Where $\alpha_{r}$ and $\alpha_{d}$ are the rise and decay temporal indices, respectively. $F_{0}$ is the normalization constant, and $t_{b}$ is the break time. The break time is related to the observed peak time ($t_{\rm p}$) by the following equation:

In the case of GRB 201015A, the smoothly broken power-law fit (smoothness parameter $k=1$) shows that initially, the optical afterglow light curve rises with a temporal index of $\alpha_{r}$=$-$1.74${}_{-0.23}^{+0.19}$ and then decays with an index of $\alpha_{d}$=1.10${}_{-0.06}^{+0.06}$, with the break time $t_{\rm b}$=217.04${}_{-18}^{+19}$ s post burst. Further, we obtained the peak time $t_{\rm p}$= 184.64${}_{-17}^{+17}$ s post burst using equation 2.

In the case of GRB 201216C, our early optical observations with FRAM-ORM telescope reveals a smooth bump in the optical afterglow light curve. We fitted the optical bump using the smoothly broken power-law function (see equation 1) and calculated the following temporal parameters: rising temporal index $\alpha_{r}$=$-0.83_{-0.41}^{+0.45}$, decay temporal index $\alpha_{d}$=$2.24_{-2.27}^{+2.14}$, and the break time $t_{\rm b}$= $161.60_{-42}^{+50}$ s post burst, respectively. Further, we determine the peak time $t_{\rm p}$= 179.90${}_{-50}^{+50}$ s post burst using the equation 2 (smoothness parameter $k$ = 3).

Reverse shock origin: According to the fireball model, the forward (moving towards the external/surrounding medium), and the reverse (propagating into the blast wave) shocks are originated due to the results of external shock between the blast-wave and ambient medium. The observed early optical peak in the afterglow light curve might be created due to the reverse shock (Kobayashi & Zhang, 2003). We used the expected temporal indices (for ISM and wind mediums in the thin and thick shells reverse shock cases) to determine the origin of early optical bump for both bursts (see Table 1 of Gao et al. (2015)). We note that the observed temporal indices values during the rising and decaying part of the bump of both the bursts are inconsistent with the expected values from the reverse shock decay in ISM or wind mediums. Therefore, the observed early bump in the optical light curves could not be due to the reverse shock.

Onset of optical afterglow: The early peak in the afterglow light curve can be produced by the onset of the afterglow and it can be used to calculate the the bulk Lorentz factor of the fireball (Sari & Piran, 1999). Ghirlanda et al. (2018) examined the early temporal coverage of GRBs with a measured redshift and calculated the bulk Lorentz factor for 67 bursts using the peak of the afterglow. In addition, they also summarize the methods used by various authors to estimate the bulk Lorentz factor. In the case of thin shell regime (T₉₀ is less than deceleration time), the peak time (in the rest frame) of the light curve provides a direct measurement of the deceleration time $t_{\rm dec}$. At the deceleration time, the Lorentz factor diminishes by a factor of two from its initial value ($\Gamma_{0}$) and enters into the self-similar deceleration phase. For the homogeneous medium surrounding the burst, the thin shell deceleration time is related to deceleration radius $R_{\rm dec}$ and bulk Lorentz factor with the following relation (Sari & Piran, 1999; Molinari et al., 2007):

Further, Ghirlanda et al. (2018) generalized the relation for both types of the surrounding medium (see equation 4).

Where s = 0 and K = 1.702 for ISM, and s = 2 and K = 1.543 for the wind-like ambient medium. $E_{\rm\gamma,iso}$ is the isotropic equivalent $\gamma$-ray energy, $m_{\rm p}$ is the proton mass, $n$ is the circumburst medium density, and $\eta$ is the radiative efficiency of the fireball.

The optical light curve of both GRB 201015A and GRB 201216C has a smooth bump that is well separated from the prompt emission, implying that the emission is in the thin shell regime. For the present analysis, we consider $\eta$ equal to 0.2 for both types of ambient medium. Further, we assume the value of $n$ = 1 cm^-3 for ISM medium. The wind medium circumburst profile is governed by the mass loss rate $\dot{M}$ and wind velocity $v_{w}$. Therefore, we used $n$ = $10^{35}\dot{M}_{-5}/v_{w,3}$ cm^-3 for wind medium (Ghirlanda et al., 2018). To determine the bulk Lorentz factor using the early bump, we have used the methodology described by (Molinari et al. 2007, hereafter M2007).

We created the SEDs using joint X-ray and optical data for both the bursts to constrain the break frequencies of the broadband synchrotron model at two epochs for GRB 201015A and one epoch for GRB 201216C, respectively. We followed the methodology given in Hu et al. (2021); Caballero-García et al. (2022); Gupta et al. (2022c) to create these SEDs. For GRB 201015A, the optical temporal index ($\alpha_{\rm o}$=$-0.92_{-0.09}^{+0.08}$) follows a simple power-law decay post optical bump and satisfies the relation $\alpha=(3p-3)/4$ in the slow cooling regime of the ISM medium, which indicates that the optical emission always remains in the $\nu_{m}$ $<$ $\nu_{o}$ $<$ $\nu_{c}$ spectral regime. The details of epochs and temporal/spectral indices are given in Table 13 of the appendix. For GRB 201015A, we performed follow-up observations with 3.6m DOT in the BVRI filters. Our 3.6m DOT optical observations covering the temporal window from 52 ks–74 ks post burst are shown in Figure 7. We created an optical SED to calculate the optical spectral index $\beta_{o2}$ and constrain the spectral regime at the epoch of DOT observations. The magnitudes are corrected for galactic extinction. In addition, we used negligible host extinction from Giarratana et al. (2022), and fitted a power-law function to determine the optical spectral index $\beta_{o2}$. We obtained a value of $\beta_{o2}$ =-1.13$\pm$0.30. The calculated value of $\beta_{o2}$ using DOT observations is shallower than ($\sim$0.3) the X-ray spectral index ($\beta_{x2}$=-1.38${}_{-0.49}^{+0.47}$) at this epoch. The measured change in the spectral indices of optical and X-ray data indicates that the cooling frequency lies in between the optical and X-ray spectral regime, although large associated errors due to limited data points and it is hard to discard other possibilities as $\beta_{o2}$ and $\beta_{x2}$ are consistent within 1 sigma. We noted that an ISM-type surrounding medium and $\nu_{m}$ $<$ $\nu_{o}$ $<$ $\nu_{c}$ $<$ $\nu_{x}$ spectral regime is also favored by the afterglow modeling of GRB 201015A by Giarratana et al. (2022) at this epoch, although, they mainly use the preliminary optical magnitudes reported in various circulars. Using the calculated value of $\beta_{o2}$, we determined the power-law index of electron distribution $p$ = (2 $\times$ $\beta_{o2}$ + 1) = 3.26 $\pm$ 0.60 at this epoch.

For GRB 201015A, we create optical-X-ray SED at two different epochs (as shown in Figure 7, upper panel). At the first epoch, the X-ray spectral index $\beta_{x1}$=-1.07${}_{-0.30}^{+0.31}$ is shallower than ($\sim$0.3) the X-ray spectral index ($\beta_{x2}$=-1.38${}_{-0.49}^{+0.47}$) at the second epoch. The evolution of X-ray spectral index $\beta_{x}$ also remain consistent within 1 sigma error. We noted that $\nu_{c}$ is just below or within the XRT band spectral regime, also favored by the afterglow modeling (Giarratana et al., 2022) at this epoch. Considering no spectral break between optical and X-ray emissions, we extrapolate the X-ray spectral index $\beta_{x1}$ toward the optical band to estimate the intrinsic optical flux. By comparing optical R-band flux with that obtained from the extrapolation, we determine an amount of extinction in the R-band $A_{R}$ = 2.2 mag. At the second epoch, the cooling frequency $\nu_{c}$ has crossed the X-ray bands at $\sim$ 52 ks.

For GRB 201216C, the late-time optical light curve seems to follow the normal decay with a power-law index, $\alpha_{o}$ = $-1.05_{-0.10}^{+0.11}$. Additionally, the X-ray temporal index $\alpha_{x}$ = $-2.21_{+0.10}^{-0.11}$ and the X-ray spectral index $\beta_{x}$ = $-0.97_{-0.05}^{+0.05}$ remain almost constant throughout the emission. The X-ray/optical temporal/spectral indices satisfied the closure relation for wind medium with emission $\nu_{m}$ $<$ $\nu_{o}$ $<$ $\nu_{c}$ $<$ $\nu_{x}$. A wind-type surrounding medium is also favored by the afterglow modeling of GRB 201216C by Rhodes et al. (2022).

For GRB 201216C, similarly, we created a optical-X-ray SED from optical to X-ray frequencies and extrapolated our X-ray spectral index (@$\sim{2.5h}$) to the lower frequency, see Figure 7, lower panel. The Galactic corrected r band magnitude lies much below the extrapolated X-ray power-law slope. Considering no spectral break and a spectral break between X-ray and optical frequencies in the SED, we estimated the intrinsic optical flux by extrapolating the X-ray spectral index at r-band frequency. By comparing the estimated optical flux to the observed VLT optical flux at 2.19 h Izzo et al. (2020), we determine the host extinction in the r-band $A_{r}\sim{8-5}$ mag, supporting the dark nature of the burst, consistent with the result of Rhodes et al. (2022).

In the present section, we compared the X-ray and the optical (see Figure 14 of the appendix) afterglow luminosity light curves of GRB 201015A and GRB 201216C with other well-studied VHE detected GRBs (GRB 160821B, GRB 180720B, GRB 190114C, GRB 190829A, and GRB 221009A). In addition, we also included X-ray (if available) and optical (if available) light curves of a nearly complete sample of nearby and supernovae-connected GRBs (GRB 050525A/SN 2005nc ($z$ = 0.606), GRB 081007A/SN 2008hw ($z$ = 0.530), GRB 091127A/SN 2009nz ($z$ = 0.490), GRB 101219B/SN 2010ma ($z$ = 0.552), GRB 111209A/SN 2011kl ($z$ = 0.677), GRB 130702A/SN 2013dx ($z$ = 0.145), GRB 130831A/SN 2013fu ($z$ = 0.479), GRB 060218/SN 2006aj ($z$ = 0.033), GRB 120422A/SN 2012bz ($z$ = 0.282), GRB 130427A/SN ($z$ = 0.340), GRB 190114C/SN ($z$ = 0.425), GRB 190829A/SN 2019oyw ($z$ = 0.078), GRB 200826A/SN ($z$ = 0.748)) for the comparison as most of VHE detected bursts are nearby and connected with supernovae. The comparison of X-ray afterglow light curves indicates that VHE detected GRB 180720B and GRB 201216C have extremely bright X-ray emission (just below the brightest X-ray emission observed from GRB 130427A at early epochs). GRB 221009A and GRB 190114C also have a very bright X-ray emission but less than those of GRB 180720B and GRB 201216C at early epochs. On the other hand, GRB 160821B, GRB 190829A, and GRB 201015A have a faint X-ray emission. The X-ray afterglow of GRB 201015A has a nearly comparable brightness with GRB 190829A. GRB 160821B, being a short burst, has the faintest X-ray light curve (after the steep decay phase) with respect to present VHE sample.

For the optical afterglow light curve comparison, we collected the R band light curve of VHE detected GRBs (other than GRB 221009A) from the literature (MAGIC Collaboration et al., 2019a; Misra et al., 2021; Gupta et al., 2021a; Fraija et al., 2019b; Jordana-Mitjans et al., 2020). For GRB 221009A, we collected the optical data from the GCN circulars¹²¹²12https://gcn.gsfc.nasa.gov/other/221009A.gcn3. For a nearly complete sample of nearby and supernovae-connected GRBs, we obtained the optical data from Kumar et al. (2022a) and references therein. We noted that GRB 180720B has the highest optical luminosity at early epochs in comparison to the present sample (similar to the X-ray light curve) and the light curve displays a smooth power-law decay across the emission period (Fraija et al., 2019b). At later phases ($\sim$ 0.2 days post burst), GRB 221009A seems to have the highest optical luminosity in comparison to the present sample. In the case of GRB 201216C, despite of the very bright X-ray emission, the optical light curve is one of the faintest with respect to present sample, typical to those observed in the case of dark GRBs. Further, we noted that GRB 201015A and GRB 201216C have a comparable optical luminous light curves. GRB 190114C (Gupta et al., 2021a), GRB 190829A (Hu et al., 2021), and GRB 201015A (Giarratana et al., 2022) had a late time bump in their optical light curve associated with their underlying supernovae explosions. In addition, we noted that GRB 201015A and GRB 201216C are the only the bursts with very early smooth bump (the onset of the forward shock) in their optical light curves with respect to present sample.

The broadband afterglow observations of the VHE detected GRBs could not be explained by the typical external shock synchrotron model (MAGIC Collaboration et al., 2019b; H.E.S.S. Collaboration et al., 2021). The multi-wavelength modeling of the observed double bump SEDs of VHE detected GRBs demands an additional SSC/Inverse Compton component. VHE detected GRBs are expected to be luminous and nearby such as GRB 180720B, GRB 190114C, and GRB 221009A. However, some of the recent detection of VHE emission from low/intermediate-luminosity bursts (such as GRB 190829A, and GRB 201015A), open a new question about their possible progenitors and viewing geometry. The high luminosity GRBs are typically observed on-axis with narrow viewing, on the other hand, low-luminosity bursts are typically observed off-axis with wide viewing angle. Recently, some authors (Sato et al., 2022; Rhodes et al., 2022) used two different jet components (narrow and wide) model to explain the origin of the observed properties of high and low luminosity VHE detected GRBs. According to this model, the early broadband emission of the high luminosity GRBs are explained using narrow jet component with typical opening angle $\theta<0.86^{\circ}$. On the other hand, the low luminosity GRBs are explained using wide jet component with typical opening angle $\theta>6^{\circ}$ (Sato et al., 2022). Sato et al. (2022) performed the broadband modeling of nearby and low-luminosity GRB 190829A using two-component jet model and noted that the prompt and afterglow emissions could not be described from the narrow jet component. The observed low/intermediate luminosity nature of GRB 190829A is explained by assuming that the viewing angle is greater than the opening angle of its narrow jet component (off-axis observations). GRB 201015A is also a nearby intermediate luminosity GRB and might have a very similar viewing geometry to GRB 190829A (viewing angle greater than the narrow jet opening angle). In the case of GRB 180720B, and GRB 190114C, the observed high luminosity nature of these GRBs are explained from emission within the narrow jet component (on-axis observations). The same applies to GRB 201216C, given its observed high luminosity.

The observed typical sub-TeV bump are explained either using SSC or external Inverse Compton (MAGIC Collaboration et al., 2019b; Abdalla et al., 2019). According to two component model, during the early phase, both the components (narrow and wide) of the jet have almost equal velocities. This negligible difference in the velocities helps to discard the possibility of the contribution from the external Inverse Compton during early phase (Sato et al., 2022). Therefore, the SSC emission mechanism is expected to explain the early VHE emission, for example GRB 190114C, and GRB 190829A MAGIC Collaboration et al. (2019b); H.E.S.S. Collaboration et al. (2021). In the case of GRB 180720B, VHE photons were observed several hours after the detection. Therefore, considering the two jets moving at different velocities, Inverse Compton significantly contributed to the observed VHE emission (Abdalla et al., 2019). In the case of GRB 201216C, considering the synchrotron process as a possible radiation mechanism at t$\sim$ 3$\times 10^{3}$ - $10^{4}$ s post burst at X-ray frequencies, we extrapolated the X-ray spectral index $\beta_{x}$ towards the GeV-TeV energies (for the spectral regime $\nu_{c}~{}<\nu_{x}$), and estimated the expected flux density $F_{\nu}\sim{2.3\times 10^{-11}}$ Jy, and $F_{\nu}\sim{3\times 10^{-13}}$ Jy at 1 GeV and 0.1 TeV, respectively. Fermi-LAT could not detect the GeV emission from GRB 201216C during the interval T₀ + 3500 s to T₀ + 5500 s post burst. This is in agreement with the fact that the expected flux density at 1 GeV at t$\sim$ 3$\times 10^{3}$ - $10^{4}$ s post burst is below the sensitivity of Fermi-LAT instrument. Furthermore, from the observed peak in the early optical light curve of GRB 201216C, we calculated the initial Lorentz factor of the burst $\Gamma\sim$ 300. With this Lorentz factor and $z=1.1$, photons of maximum energy, $\leq$15 GeV, are allowed by synchrotron process Fraija et al. (2019b).

Similarly, for GRB 201015A expected flux density $F_{\nu}\sim{3.8\times 10^{-13}}$ Jy, and $F_{\nu}\sim{2.8\times 10^{-15}}$ Jy at 1 GeV and 0.1 TeV, respectively, is obtained at 3ks-5ks. A bulk Lorentz factor $\sim$200 is obtained from the observed peak in the optical light curve. With this Lorentz factor and $z=0.426$, we estimated the maximum energy of synchrotron photons $\leq$14 GeV Fraija et al. (2019b).

Therefore, observed early VHE emission could not be explained using the synchrotron emission model. The previous discussion suggests that the early detection of VHE photons from GRB 201015A and GRB 201216C required that the ultra-relativistic outflow boosts the energy of low-energy synchrotron photons to the VHE energies via the SSC process.