GRB 210121A: Observation of photospheric emissions from different regimes and the evolution of the outflow

GRB 210121A was observed by Insight-HXMT, Gravitational wave high-energy Electromagnetic Counterpart All-sky Monitor (GECAM), Fermi Gamma-ray Burst Monitor (Fermi/GBM) on Jan 21st 2021. In this work, photospheric emission from a structured jet is preferred to interpret the prompt emission phase of GRB 210121A and emissions from different regimes are observed on-axis. Particularly, the emission from the intermediate photosphere is first observed in the first 1.3 s of the prompt emission, while those from the other part are dominant by the saturated regime, and offers an alternative explanation compared with the previous work. Moreover, the emissions with considering the intermediate photosphere can well interpret the changes on low-energy photon index $\alpha$ during the pulses. Besides, the evolution of the outflow is extracted from time-resolved analysis, and a correlation of $\Gamma_0 \propto L^{0.25\pm0.05}_0$ is obtained, which implies that the jet may be mainly launched by neutrino annihilation in a hyper-accretion disk.


INTRODUCTION
Despite about a half century of observations, prompt gamma-ray burst (GRB) emission mechanisms are still a matter of interest. Two leading scenarios have been suggested to interpret the observed spectra of GRBs. One is synchrotron radiation, which invokes a non-thermal emission of relativistic charged particles either from internal shocks or from internal magnetic dissipation processes (Lloyd & Petrosian 2000;Tavani et al. 2000; Baring & Braby 2004;Zhang 2018). The fast cooling problem in relevance to the lowenergy photon index α could be relieved by synchrotron radiation of electrons in a moderately fast-cooling regime (Uhm & Zhang 2014) or detailed treatment of the cooling electrons (Derishev et al. 2001;Geng et al. 2018). As a natural consequence of the fireball model, photospheric emission produced from highly relativistic outflows was previously considered as an explanation for prompt gammaray bursts (Goodman 1986;Paczynski 1986), where the optical depth at the base of the outflow is much larger than unity (Piran 1999). The Planck spectrum related to the photospheric emission could be broadened in two ways. First, dissipation below the photosphere can heat electrons above the equilibrium temperature. These electrons emit synchrotron emission and comptonize the thermal photons, thereby modify the shape of Planck spectrum (Pe'er et al. 2005;Pe'er et al. 2006;. Observational evidence for subphotospheric heating has been provided by Ryde et al. (2011). Besides, internal shocks below the photosphere , magnetic reconnection (Thompson 1994;Giannios & Spruit 2005), and hadronic collision shocks (Beloborodov 2010;Vurm et al. 2011) can also cause dissipation. Second, the modification of Planck spectrum could be caused by geometrical broadening. Photospheric radius is found to be a function of the angle to the line of sight of the photons of thermal emission observed (Abramowicz et al. 1991;Pe'er 2008). This means that the observed spectrum is a superposition of a series of blackbodies of different temperature, arising from different angles to the line of sight. Moreover, Pe'er (2008) showed that photons make their last scatterings at a distribution of radii and angles. The observer sees simultaneously photons emitted from a large range of radii and angles. Therefore, the observed spec-trum is a superposition of comoving spectra (Ryde et al. 2010;Hou et al. 2018). Lundman et al. (2012) studied the non-dissipative photospheric (NDP) emissions from a structured jet, and reproduced the average low-energy photon index (α = −1) independent of viewing angle. The observed evolution patterns of the νF ν peak energy (E p ), including hard-to-soft and intensity-tracking could be reproduced by this model as well (e.g. Deng & Zhang 2014;Meng et al. 2019).
There are three regimes to be discussed in photospheric emission from a structured jet: (I) unsaturated emissions which are dominant by the regime of unsaturated acceleration (the saturation radius is greater than the photonsphere radius, R s > R ph ); (II) saturated emissions from the regime of R s ≤ R ph works for all over the wind profile; (III) the intermediate photosphere (Song & Meng 2022) which represents the case where the regimes of R s > R ph and R s ≤ R ph work in lower and higher latitudes respectively, and the contribution from the latter can not be ignored. It has been never mentioned in the previous study in the spectral fit to photospheric emissions. Besides of the off-axis NDP model in unsaturated regime (Wang et al. 2021), in this paper, we find that the on-axis NDP model with considering intermediate photosphere is an alternative description for the emissions of GRB 210121A.
The evolution of the outflow is also extracted, to offer an interpretation about the mechanism for launching the jet. If we take the hyper-accreting black hole (BH) as the central engine, the GRB jet may be launched from two mechanisms. One is νν annihilation in a neutrinodominated accretion flow (NDAF) (e.g. Popham et al. 1999;Narayan et al. 2001;Matteo et al. 2002;Kohri & Mineshige 2002;Gu et al. 2006;Chen & Beloborodov 2007;Janiuk et al. 2007;Lei et al. 2009;Liu et al. 2010), and generates a fireball which is dominant by the thermal component. The jet is launched by neutrino annihilation (νν → e + e − ), andĖ νν is the neutrino annihilation power. The other one is Blandford&Znajek mechanism (BZ, Blandford & Znajek 1977). The spin energy of the BH is tapped by a magnetic field, and produces a Poynting flux. The correlations between the baryon loading parameter and the power for these two mechanisms are both positive, η ∝Ė 0.26 νν (Lü et al. 2012) for the former, and µ 0 ∝Ė 0.17 BZ (Yi et al. 2017) (µ 0 = η(1 + σ 0 ), where σ 0 is the ratio of Poynting flux luminosity to the matter flux) for the latter. Thus, the index of the correlation offers a criterion.
This paper is organized as follows. In Section 2, the observation of GRB 210121A by different missions is introduced. In Section 3, the NDP model and R ph in different regimes are introduced, especially in the in-termediate photosphere. In Section 4, the methods for binning, background estimation and the spectral fitting are clarified. In Section 5, time-resolved analyses are performed; the properties of the jet are extracted, while the regimes are determined. Then the discussion and conclusion are given in Section 6. The former one is taken to be the T 0 in the following analysis. The photon flux of GRB 210121A is shown in Figure 1, which is extracted from the joint analysis utilized in . The cutoff power law (CPL) model is preferred for almost all the time-resolved slices as shown in Wang et al. (2021). The values of α of the first epoch obtained from time-resolved analysis are greater than -2/3, as well as those of the most part of the second epoch from T 0 +2.8s to T 0 +14.9s. According to Meng et al. (2021), the NDP spectra from a hybrid outflow with a moderate magnetization have a larger high-energy photon index (β), and are more compatible with the BAND function rather than the CPL model. Therefore, we prefers a pure hot fireball (or the Poynting flux completely thermalized below the photosphere) in the prompt phase.
We can not justify the radiation mechanism only from α (Burgess et al. 2020;Meng et al. 2018), thus a physical model fitting is needed.

THE MODELING OF PHOTOSPHERIC EMISSIONS FROM DIFFERENT REGIMES
A structured jet is supported by observation, (e.g., studies on GRB 170817A, Lazzati et al. 2017;Bromberg et al. 2018;Geng et al. 2019) and the relativistic magnetohydrodynamic simulations for the GRB jet are also performed (e.g., Kathirgamaraju et al. 2019). In this analysis, we assume a structured jet with an opening angle θ c and luminosity L 0 . Motivated by the results of MacFadyen & Woosley (1999) and Zhang et al. (2003) e.g., the jet is structured with an innerconstant and outer-decreasing angular baryon loading parameter profile with the form for the GRB prompt emission phase (Dai & Gou 2001;Rossi et al. 2002;Zhang & Mészáros 2002;Ku- where η is the angle-dependent baryon loading parameter which is also the bulk Lorentz factor Γ in the saturated acceleration regime; η 0 is the maximum η, and also denoted as Γ 0 ; θ is the angle measured from the jet axis; θ c is the half-opening angle for the jet core; p is the power-law index of the profile; η min = 1.2 is the minimum value of η, differing from unity for numerical reasons. The exact value of η min only affects the very low energy spectrum, many orders of magnitude below the observed peak energy (Lundman et al. 2012). This baryon loading parameter profile have used in Lundman et al. (2012), Meng et al. (2018) and Meng et al. (2019). Flux of observed energy E obs at the observer time t in the case of the continuous wind is deduced and explained in Section A.1. The angle-dependent photosphere radius R ph , as the radius from which the optical depth for a photon that propagates in the radial direction is equal to unity, is defined as where β is the velocity, r 0 is the radius of the central engine and dṀ (θ)/dΩ = L 0 /4πc 2 η(θ) is the angle-dependent mass outflow rate per solid angle; m p is the mass of the proton, c is the light speed, and σ T is electron Thomson cross section. Note the unsaturated regime used in Wang et al. (2021) is R ph R s . In this paper, an unsaturated emissions of R ph R s is considered and added between the case of R ph R s and R ph ≥ R s (Song & Meng 2022).
R s = η(θ)r 0 monotonically decreases with θ, while R ph monotonically increases with θ. Thus, there exists a critical value θ cri , and for θ ≥ θ cri , it satisfies R s (θ) ≤ R ph (θ), namely, Note that for the intermediate photosphere, θ cri satisfies 0 < θ cri < 5/Γ 0 , and the contribution from the saturated regime contributes to the prompt emission (Lundman et al. 2012), or the unsaturated emission is dominant.
For an intermediate photosphere, R ph is described by two forms: R ph in lower latitude of the unsaturated part is described by the second item of Equation (2), while that of higher latitude takes the saturated form.
If for a set of parameters, θ cri is found to be 0 by Equation (3), the saturated regime works for all over the jet profile. Therefore in our analysis, the regime could be determined by the parameters obtained from the fitting to spectra. In this work, the Bayesian blocks (BBlocks) method introduced by Scargle et al. (2013) and suggested by Burgess (2014), is applied with a false alarm probability p 0 = 0.01 on light curves. In some cases, the blocks are coarse for fine time-resolved analysis. Burgess (2014) suggested that the constant cadence (CC) method is accurate when the cadence is not too coarse. Therefore, we take a combination of BBlocks and CC methods, fine binning of constant cadence are performed in each block, and only the bins with the signal-to-noise ratio (S/N) ≥ 20 at least in one detector should be utilized. 15 bins from [T 0 -0.01, T 0 + 14.90] s are shown and labeled in Figure 1, where 5 bins are in the first epoch and the others are in the second one.

Background estimation and spectral fitting method
A polynomial is applied to fit all the energy channels and then interpolated into the signal interval to yield the background photon count estimate for GRB data. Markov Chain Monte Carlo (MCMC) fitting is performed to find the parameters with maximum Poisson likelihood. The Sampling tool is emcee (Foreman-Mackey et al. 2013). Wei et al. (2016) suggested a bayesian information criterion (BIC) as a tool for model selection, a model that has a lower BIC value than the other is preferred. If the change of BIC between these two models, ∆BIC is from 2 to 6, the preference for the model with the lower BIC is positive; if ∆BIC is from 6 to 10, the preference for that is strong; and if ∆BIC is above 10, the preference is very strong.

FIT RESULTS
In this section, time-averaged spectral fitting for two epochs as well as time-resolved analyses for fine bins are performed with the NDP model. The properties of the jet and the evolution of the outflow are extracted. 5.1. time-averaged results of two epochs Table 1 and Figure 2 show the fit results, MCMC samples and spectra for the two epochs. It is found that θ c is at 0.01 with small uncertainties for both epochs. The values of r 0 are well consistent with 10 7.61 cm. p becomes larger in epoch 2 than that in epoch 1. The luminosity of outflow changes with time, and becomes smaller in epoch 2 than that in epoch 1, which almost has a similar trend to the flux. For the first epoch, χ 2 =248.2 and BIC=262.14 with the degree of freedom (dof) of 203, while χ 2 =360.0 and BIC=375.15 with the dof of 323 in the second epoch. We also performed the fit with an on-axis NDP model dominant by the regime of R ph R s . However, the fit is not good, with BIC equals to 904.85 and 2014.26 for the two epochs with the same degrees of freedom. Thus, for on-axis observation, the NDP model in the regime of R ph R s is not preferred, and inconsistent with the spectra.
These two results corresponds to θ cri = 0, which means the emission from the saturated regime are dominant. If considering with changes on L 0 and η 0 during the epochs, the emissions may be from the intermediate photosphere, or the unsaturated regime, which are determined in the fine time-resolved results.

time-resolved results for fine bins
Since the uncertainties of fit results from timeaveraged spectra seem not large, we could assume that r 0 , p, θ c have small changes during each epoch; z stays constant for the whole GRB duration. Therefore, they could be fixed in the time-resolved fitting to suppress the uncertainty in extracting the correlation of Γ 0 − L 0 . The timeresolved analysis is performed with float η 0 and L 0 , while other parameters are fixed to p = 1.0 for bins 0-4 in the first epoch and p = 1.27 for bins 5-14 for the second epoch. For both epochs, logr 0 =7.61, θ c = 0.01, z=0.37. The fit results are listed in Table 1, where four bins have non-zero θ cri . All over the bins, bin 1 has the largest η 0 and flux, as well as the largest θ cri . We can conclude that it has larger contribution from the unsaturated regime than the other bins. The relation of θ cri and L 0 for each set of (r 0 , p, θ c ) of bin 1 and bin 5 are shown in Figure 3 in red dotted lines. The values of luminosity measured in time-resolved analysis are represented by vertical dot-dashed lines in blue. As shown in Figure 3 (a), the corresponding θ cri of bin 1 denoted as a red star is less than θ c (denoted by the green solid line) as well as 5/Γ 0 (denoted by the orange dashed line), which means that the emission is from the intermediate photosphere. For comparison, emissions in other bins, e.g., bin 5, is shown in Figure 3 (b). The blue vertical line does not intersect with the red dotted line representing θ cri , which means the corresponding θ cri is 0, and the emission is from the saturated regime. The time-resolved spectra are shown in Figure 4 (a), where the luminosity and hardness decrease by time generally.
The relations between η 0 and L 0 are shown in Figure 4. Γ 0 ∝ L 0.25±0.05 0 is extracted from the timeresolved results denoted by black dots. The timeaveraged results for two epochs represented by red squares are also plotted, and consistent well with this trend.

DISCUSSION AND CONCLUSION
In the previous study, Wang et al. (2021) performed a fit on the first epoch in GRB 210121A with the NDP model in unsaturated regime (R ph R s ), and the observation is determined to be greatly off-axis with the extracted θ v 15/Γ 0 . In this paper, we use NDP model with considering the intermediate photospheric emissions to describe two epochs of GRB 210121A and give an alternative description with an on-axis observation. The obtained luminosity, log(L 0 (erg s −1 )) = 50.62 +0.42 −0.26 , is naturally smaller than that of off-axis, log(L 0 (erg s −1 )) = 51.94 +0.84 The intermediate photosphere always has a moderate α between these two regimes (Song & Meng 2022), where α is extracted from spectra below the peak energy with an exponential cut-off power law. This could explain that bin 1 has the largest α = −0.22 +0.07 −0.07 (this value is from Wang et al. (2021)), while the others has smaller θ cri and more contribution from the saturated regimes. According to Wang et al. (2021), it seems that the off-axis unsaturated NDP model could also give a smaller α than on-axis unsaturated NDP model, therefore, there may be double solutions for the observation, and the on-axis NDP model with considering intermediate photonsphere gives an alternative solution. Liang et al. (2010) presented a correlation of Γ 0 ∝ E 0.25 iso,γ . Lü et al. (2012) discovered an even tighter correlation Γ 0 ∝ L 0.30 iso,γ from 50 GRBs and proposed an interpretation. Considering the beaming factor (f b ∝ L −0.145 iso,γ ) and the similar efficiency of jet kinetic energy via internal shock in dissipation to prompt emissions of GRBs, it has Γ 0 ∝ L 0.22 iso,γ . However, Yi et al. (2017) found that the data are more consistent with the latter mechanism, the index∼ 0.14. In this analysis, Γ 0 −L 0 of GRB210121A is obtained from the time-resolved anal- ysis, and there are some advantages. First, we do not need to consider the efficiency of the prompt emission, and L 0 ∝Ė (Ė is the power). Second, we obtain the baryon loading parameter directly from the the fit result, and consider the regime in this procedure. The Lorentz factor of the outflow may be lower than the baryon loading parameter if the emission is from the unsaturated regime. η 0 ∝ L 0.25±0.05 0 is extracted from time-resolved analysis, and it is more consistent with the mechanism of νν annihilation.
In summary, the NDP model with considering the intermediate photosphere is first extracted in GRB 210121A, and it could well explain the changes on α during the pulses. The correlation of η 0 −Γ 0 is extracted from the time-resolved analysis, implies the jet of GRB 210121A may be mainly from the νν annihilation.   As discussed in (Lundman et al. 2012;Deng & Zhang 2014;Meng et al. 2018, e.g.), flux of observed energy E obs at the observer time t in the case of continuous wind could be shown as in Equation (A1), where the velocity β = v c and the Doppler factor D = [Γ(1 − β cos θ LOS )] −1 both depend on the angle θ to the jet axis of symmetry, in which θ LOS is the angle to the line of sight (LOS) of the observer. The viewing angle θ v is the angle of the jet axis of symmetry to the LOS. If it is on-axis, we have θ v = 0 and θ LOS = θ. dṄ γ /dΩ =Ṅ γ /4π anḋ N γ = L/2.7k B T 0 , where L is the total outflow luminosity, T 0 = (L/4πr 2 0 ac) 1/4 is the base outflow temperature and a is the radiation constant. dP/dE describes the probability for a photon to have an observer frame energy between E and E + dE within volume element dV , and it is a comoving Planck distribution with the comoving temperature T (r, Ω), as shown in Equation (A2), where k B is the Boltzmann constant, T ob (r, Ω) = D(Ω)· T (r, Ω) is the observer frame temperature, and in the saturated acceleration regime ( R s < R ph ), it is defined by ( , R s (Ω) < R ph (Ω) < r.

(A3)
In the unsaturated acceleration case, comoving temperature T (r) is given by where Γ is R ph /r 0 . The spectra are always obtained within a time interval in time-averaged or time-resolved analysis. The luminosity is taken to be constant as an averaged value, and the spectra is taken as a time-integrated spectra in this time interval. Continuous wind could be assumed to consist of many thin layers from an impulsive injection at its injection timet, and the wind luminosity att is denoted as L w (t). The flux of observed energy E obs at the observer time t is given by where t D is the duration of emission of the central engine, and r = βc(t−t) u in F obs E (θ v , t,t, L w (t)) in Equation (A1). The observed GRB has a duration, thus, it is reasonable that the central engine produces a continuous wind. For the case of the constant wind luminosity, after the central engine has an abrupt shutdown (t > t D ), the flux sharply drops. Moreover, the observed light curves of the GRBs show relatively slow change in luminosity, unlike the steep rise and fall (almost within 10 −3 s), thus, a continuous wind with variable luminosity is used to simulate a GRB pulse.