{'fraction_non_alphanumeric': 0.04083226418251799, 'fraction_numerical': 0.03190717911530094, 'mean_word_length': 4.167867781301047, 'pattern_counts': {'":': 0, '<': 23, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Low-energy Spectra of Gamma-ray Bursts from Cooling Electrons

The low-energy spectra of gamma-ray bursts' (GRBs) prompt emission are closely related to the energy distribution of electrons, which is further regulated by their cooling processes. We develop a numerical code to calculate the evolution of the electron distribution with given initial parameters, in which three cooling processes (i.e., adiabatic, synchrotron and inverse Compton cooling) and the effect of decaying magnetic field are coherently considered. A sequence of results are presented by exploring the plausible parameter space for both the fireball and the Poynting-flux-dominated regime. Different cooling patterns for the electrons can be identified and they are featured by a specific dominant cooling mechanism. Our results show that the hardening of the low-energy spectra can be attributed to the dominance of synchrotron self-Compton cooling within the internal shock model, or to decaying synchrotron cooling within the Poynting-flux-dominated jet scenario. These two mechanisms can be distinguished by observing the hard low-energy spectra of isolated short pulses in some GRBs. The dominance of adiabatic cooling can also lead to hard low-energy spectra when the ejecta is moving at an extreme relativistic speed. The information from the time-resolved low-energy spectra can help to probe the physical characteristics of the GRB ejecta via our numerical results.


INTRODUCTION
The radiation mechanism responsible for the prompt emission of gamma-ray bursts (GRBs) remains unidentified since the discovery of GRBs. A typical spectrum of the GRB prompt emission can usually be well fit by the so-called Band function (Band et al. 1993), which smoothly joins low-and high-energy power laws. Except for the Band component, the superposition of multiple spectral components was also observed in some GRBs, such as the thermal component (e.g., Ghirlanda et al. 2003;Ryde 2005;Ghirlanda et al. 2007;Ryde et al. 2010;, or an additional power-law component (e.g., Abdo et al. 2009a;Ackermann et al. 2011). Since the Band function is still an empirical description for the GRB spectra, further studies are needed to manifest its physical origin.
Synchrotron radiation of electrons has been suggested to be a possible mechanism. However, one problem (called fast cooling problem, see Ghisellini et al. 2000;Zhang & Yan 2011) remains in a simple synchrotron model, i.e., the low-energy spectral index α (F ν ∝ ν α ) is predicted to be −1/2 for fast cooling electrons (Sari et al. 1998), which is incompatible with the fact that the observed value is ∼ 0 in the majority of GRBs (Band et al. 1993;Preece et al. 2000;Nava et al. 2011;Geng & Huang 2013). Modified synchrotron models have been proposed to ease this conflict. When electrons cool mainly via inverse Compton (IC) scattering in the Klein-Nishina (KN) regime, it was suggested that the energy loss rate of electrons is roughly ∝ γ −1 e (γ e is the electron Lorentz factor), and the corresponding flux density is F ν ∼ ν 0 (see Derishev et al. 2001;Bošnjak et al. 2009;Wang et al. 2009;Nakar et al. 2009;Fan 2010;Daigne et al. 2011). Subsequent detailed analytical study on this solution shows that it is impossible to obtain a spectrum with α > −0.1 using IC cooling in the KN regime (Barniol Duran et al. 2012). Recently, by solving the continuity equation of electrons in energy space numerically, Uhm & Zhang (2014) found that the fast cooling electrons could have a harder energy spectrum when the surrounding magnetic field is decreasing. Their numerical results show α can be even harder than 0 in some certain parameter regime. However, IC cooling of electrons was not included in their calculations. It is then crucial to take IC cooling into account when solving the continuity equation of electrons.
Alternative models based on photospheric emission have also been proposed to explain the GRB prompt emission (e.g., Rees & Mészáros 2005;Giannios & Spruit 2007;Beloborodov 2010;Ryde et al. 2011). Indeed, the spectra of some GRBs are found to be consistent with a photospheric component (e.g., Ryde et al. 2010;Larsson et al. 2011;Pe'er et al. 2012). However, the photosphere model typically predicts α ∼ 1.4 (Deng & Zhang 2014). Some additional effects should be considered to explain the observed index of α ∼ 0 with the photosphere model (e.g., Lundman et al. 2013). In general, given that the main spectral component of a arXiv:1709.05899v3 [astro-ph.HE] 30 Nov 2017 typical burst is the Band component, it is still rational to suppose that the emission comes from a non-thermal mechanism in an optically thin region Veres et al. 2012;Kumar & Zhang 2015). So in this study, we work within the framework of synchrotron radiation and focus on the low-energy spectra of GRBs.
In principle, the profile of synchrotron spectra of GRBs observed is directly determined by the distribution of electrons (or called the electron spectrum/distribution for short) in the ejecta. The initial spectrum of the electrons that are accelerated somehow will soon be modified by cooling processes, thus electron cooling is a key factor in the prompt emission, especially when we focus on the low-energy spectra of GRBs. Although only synchrotron emission is observed by us in the keV -MeV band, the electrons actually can be cooled in three ways, which are adiabatic, synchrotron, and IC cooling respectively. Numerical solution of the electron distribution considering the three processes has been presented in previous researches within the internal shock scenario (Bošnjak et al. 2009;Daigne et al. 2011). They found the majority of observed GRB prompt spectra can be reconciled with a synchrotron origin. Some useful constraints on the microphysics of internal shocks were presented. This numerical approach can also be extended to the Poyntingflux-dominated jet scenario. In this study, we investigate electron cooling in different physical situations and give a clue to distinguish them.
In addition to modeling the time-integrated spectra analytically, the analyses on the time-resolved spectra of the prompt emission (e.g., Lu et al. 2012;Zhang et al. 2016;Jiang et al. 2016) can also provide important information on the radiation process. On the other hand, the numerical method has the advantage over the analytical method in that it can incorporate different radiation mechanisms and can follow the evolution of physical properties in the emitting region. For example, Daigne et al. (2011) predicted that the high-energy (> 100 MeV) lightcurve may display a prolonged pulse duration due to the IC emission. Uhm & Zhang (2014) revealed that as a jet expands rapidly from the central engine, the magnetic field in the emission region decreases, resulting in harder (than the case of constant magnetic field) spectra. Therefore, in our study, we develop a numerical code to calculate the evolution of the electron spectra and the corresponding flux spectra with different parameter sets. Adiabatic, synchrotron, IC cooling of electrons (also see Bošnjak et al. 2009) and the geometric effect of the emitting shells are considered properly in our code. Using this code, we can explore the resulting spectra in the plausible parameter space, which may provide clues to help relate the observed GRB spectra with the physical processes in the GRB ejecta. Different kinds of cooling patterns for electrons obtained in different scenarios can also serve as the baseline for further explorations.
The structure of this article is as follows. The three main cooling processes considered are briefly described in Section 2. The constraints from observations on the parameters involved in our calculations are presented in Sections 3. In Section 4, we derive the conditions under which one particular cooling process will be dominant analytically. The analytical results are then compared with the numerical results in Section 5, where we generally study the roles played by different cooling processes in determining the evolution of the low-energy electron distribution. Finally, in Section 6, we summarize and discuss our results. The details of the numerical method and some relevant formulations used are given in Appendices A and B respectively.

COOLING OF ELECTRONS
In the co-moving frame of a relativistic jet, when an electron with a Lorentz factor of γ e is moving in the magnetic field of strength B , it will lose energy by the synchrotron radiation at a rate of (Rybicki & Lightman 1979) where σ T , m e , and c are the Thomson cross-section, electron mass, and speed of light respectively. Hereafter, the superscript prime ( ) is used to denote the quantities in the co-moving frame. The electron also undergoes adiabatic cooling (Uhm et al. 2012;Geng et al. 2014), i.e., where we have taken the co-moving electron number density n e ∝ R −2 for an expanding shell. Additionally, the electrons will be cooled by the IC scattering of self-emitted synchrotron photons, which is referred to as synchrotron self-Compton (SSC) process. The SSC cooling rate is given by (Blumenthal & Gould 1970;Fan et al. 2008) γ e,SSC = − 1 m e c 2 3σ T c 4γ 2 e ν max ν min n ν dν ν ν ic,max ν ic,min hν ic dν ic F (q, g), (3) where F (q, g) = 2q ln q + (1 + 2q)(1 − q) + 1 2 (4qg) 2 1+4qg (1 − q), g = γ e hν mec 2 , w = hν ic γ e mec 2 , q = w 4g(1−w) . The upper limit of the internal integral can be derived as hν ic,max = γ e m e c 2 4g 4g+1 , and the lower limit is ν ic,min = ν . Overall, the total cooling rate of an electron can be obtained by summing up the processes mentioned above, i.e., γ e,tot =γ e,syn +γ e,adi +γ e,SSC .
Heating of low-energy electrons due to synchrotron absorption (Ghisellini & Svensson 1991;Gao et al. 2013) is not considered here, which may pile up electrons in the low-energy range. In this work, we focus on the cooling processes in order to investigate their roles clearly. Heating or acceleration of electrons will be incorporated in our future studies. The GRB prompt emission comes from a group of electrons, of which the instantaneous spectrum can be denoted as dN e /dγ e . This electron distribution can be obtained by solving the continuity equation of electrons in energy space (Longair 2011) ∂ ∂t dN e dγ e + ∂ ∂γ e γ e,tot dN e dγ e = Q(γ e , t ), where Q(γ e , t ) is the source function that describes the electrons injected into the emitting region. If the bulk Lorentz factor of the jet is Γ, the co-moving time t can be related to the observer's time by and the position of the jet head is described by where β = 1 − 1 Γ 2 is the dimensionless velocity of the jet.

CONSTRAINTS FROM OBSERVATIONS
On one hand, the characteristics of the emission site for the prompt emission, the relativistic jet (or the relativistic ejecta), are still under research. On the other hand, the knowledge on the composition of the jet, and some quantities such as the bulk Lorentz factor of the jet, the emission radius, the strength of B , and the Lorentz factor of the electrons are crucial to model the GRB spectra. Here, rather than assuming a detailed model, we try to derive some constraints on relevant quantities from observations. These constraints have been derived analytically in prior articles (e.g., Kumar & McMahon 2008;Beniamini & Piran 2013Kumar & Zhang 2015). The main logic of our derivation is similar to these works. With proper ranges of these quantities, we can then analyze the cooling behavior of electrons and correspondingly perform calculations in the following sections.
Assuming that the bulk Lorentz factor of the GRB ejecta is Γ, the Lorentz factor of electrons that radiate at the GRB spectral peak energy E peak is γ m , then we have where h, q e are the Planck constant and electron charge respectively, z is the redshift of the burst. The radiative cooling time for an electron of γ e in the observer frame is where Y is the Compton-Y parameter. The dynamical time of the jet can be expressed as t d ∼ R(1 + z)/(2Γ 2 c). For typical parameters, the magnetic field strength in the emission region is strong enough that the electrons are in the fast cooling regime (also see Equation (13)). The fast cooling condition requires that t c (γ m ) ≤ t d , which implies Taking typical values of E peak 500 keV, and R 10 15 cm, we can get, B Γγ 2 m = 2.9 × 10 13 (1 + z) Γγ m ≤ 3.3 × 10 7 (1 + z) 2/3 (1 + Y ) 1/3 E peak 500 keV 2/3 R 10 15 cm 1/3 , by using Equations (8) and (10).
On the other hand, the specific flux at E peak in the observer frame can be expressed as (e.g., Beniamini & Piran 2013Kumar & Crumley 2015) where N e is the total (already corrected for 4π solid angle) number of electrons with γ e > γ m , and D L is the luminosity distance of the burst. Then we can estimate the number of electrons needed to produce a given observed flux by combining Equations (11) and (13), The corresponding average injection rate of electrons is where we have used t c (γ m ) as the injection timescale in order to maintain the electron distribution and the intensity of the radiation flux within this period. One can see that in a synchrotron model, the typical value of N inj is almost model independent and may be compared/verified with further detailed simulation results. Now, we consider two leading models respectively to obtain the plausible range of the key parameter γ m . First, if the jet is magnetically dominated, i.e., a relativistic Poynting-flux-dominated jet 9 (e.g., Zhang & Yan 2011;Beniamini & Piran 2014;Kumar & Crumley 2015), its isotropic equivalent magnetic luminosity is L B B 2 8π Γ 2 4πR 2 c. The kinetic energy power of the accelerated electrons is L e N inj m e c 2 γ m Γ 2 . In principle, the ratio of L e to L B , should be less than 1 since the jet is magnetically dominated. The upper limit of γ m can thus be derived from Equation (16). However, one should note that η e may still be slightly larger than 1 in the realistic case, since the magnetic field B that cools electrons here may be smaller than the average magnetic field of the ejecta.
In the framework of the internal shock model (e.g., Rees & Meszaros 1994;Daigne et al. 2011), the dominant energy of the jet should be the kinetic energy carried by protons. If we consider that there are η p protons 9 In this article, by saying the scenario of Poynting-fluxdominated jet, we mean the regime invoking a large emission radius (the magnetization parameter is not necessarily very large, see Zhang & Yan 2011), which is consistent with a Poynting-fluxdominated regime. Note that besides the synchrotron mechanism, emission from Poynting-flux-dominated jet has also been discussed in the photosphere context (e.g., Drenkhahn & Spruit 2002;Metzger et al. 2011). for every accelerated electron and assume that the accelerated protons remain non-relativistic, then we get the kinetic energy power of protons as L p η p N inj m p c 2 Γ 2 (Bošnjak et al. 2009;Daigne et al. 2011;Beniamini & Piran 2013), and the ratio between L B and L p is Here, ξ B can be equivalently treated as the familiar magnetization parameter σ at the emission radius R. On the other hand, the ratio between L e and L p is This gives an upper limit of γ m , i.e., γ m = ξ e ηpmp me < ηpmp me (Barniol Duran et al. 2012;Beniamini & Piran 2013;Kumar & Crumley 2015) since ξ e should be less than 1. Combining this limit with Equation (17), we obtain the plausible range of ξ B , Consequently, the range of γ m could be derived from Equations (17) and (19). One may notice that the free parameter η p is crucial to determine the ranges of other quantities in the internal shock model. Previous studies indicate that η p ≥ 10 should be satisfied for the internal shock model to explain the GRB spectra (see Kumar & Zhang 2015 for a review). Whether this requirement can be fulfilled within the simulation of collisionless ionelectron shocks is still under debate (Sironi et al. 2015). This issue goes beyond the scope of our current study. In this work, we admit η p ≥ 10 first and see whether the low-energy spectra can be explained naturally.

DIFFERENT REGIMES
With the estimates on the ranges of the key parameters shown above, we now discuss the possible cooling behaviors of electrons in the emission region of a GRB analytically. Conventionally, we first take the synchrotron radiation as the main cooling process for electrons, since the spectra of the observed prompt emission resemble the synchrotron spectra. However, it is possible that the electrons could also lose energy largely through the other two processes. For instance, the adiabatic cooling rate for an electron of γ e will dominate the synchrotron cooling rate ifγ e,adi γ e,syn ≥ 1, which further gives 5.3 × 10 22 (1 + z) 2 E peak 500 keV 2 R 10 15 cm by using Equations (1), (2) and (11). For electrons of γ e ≤ 10 3 , one would find this situation occurs when Γ ≥ 10 3 and γ m ≥ 10 4 , or when R is significantly smaller than 10 15 cm. We will see that adiabatic cooling does dominate in some cases in the following calculations.
The corresponding physical requirements are not straightforward sinceγ e,SSC involves double integral in Equation (3). Before giving accurate numerical results, some simple estimates can be done primarily. The magnetic energy density in the co-moving frame of the ejecta is while the radiation energy density in the co-moving frame can be calculated as If scattering between an electron of γ e with photons is always in the Thomson regime, it is well known thaṫ γ e,SSC /γ e,syn can be approximated as U γ /U B . We define γ T as the Lorentz factor of the electrons above which the scattering with the E peak photons is in the KN regime (Wang et al. 2009;Nakar et al. 2009), i.e., γ T will be smaller than γ m only if Γ is not very large. So, for electrons within the range of γ T < γ e < γ m , IC cooling is in the KN regime.
We also define hν KN = m e c 2 /γ e as the critical photon energy of which the IC scattering between the electron of γ e is in the KN regime. The scattering between the electron with photons of frequency ν < ν KN can effectively cool the electron. Assuming a low-energy photon spectrum of νF ν ∝ ν δ (hν < E peak ), then we havė  (14) into Equation (26). If the system is steady (∂/∂t = 0), the lower limit of δ is 0.5 when electrons are in fast cooling due to the synchrotron radiation (dN e /dγ e ∝ γ −2 e for γ e < γ m ), while the upper limit is 1 when electrons are cooled by the SSC radiation (dN e /dγ e ∝ γ −1 e for γ e < γ m , Wang et al. 2009). One will then find that for γ m ≥ 10 4 , the condition in Equation (27) can be met at least for electrons of γ e γ T . Moreover, for even smaller emission radius (e.g., R ≈ 10 14 cm), this condition will be relaxed significantly. So it is essential to consider SSC cooling during the evolution of the electron distribution. In summary, Equations (21) and (27) obtained are useful explicit criteria on judging how an electron cools with given relevant parameters.

NUMERICAL CALCULATIONS
The main task is to solve the continuity equation of electrons in the energy space, i.e., Equation (5), which is also called the advection equation with a source term. This kind of partial differential equations can be efficiently solved by the constrained interpolation profile (CIP) method (Yabe & Aoki 1991;Yabe et al. 2001). Detailed discretization procedure can be found in Appendix A. In principle, the final results are determined by the initial and boundary conditions for Equation (5). On the other hand, we have already obtained the plausible range for relevant parameters according to the estimates in Section 3. As we mentioned before, we intend to give an overview of the evolution of the electron distribution under different physical conditions. So we explore the parameter space by performing several groups of calculations to investigate the roles played by different radiation mechanisms in different cases. In this paper, we adopt the assumption that the co-moving magnetic field in the jet is decaying with radius as proposed in (Uhm & Zhang 2014), i.e., where B 0 is the magnetic strength at R 0 , and R 0 is the radius where the jet begins to produce the first photon that observed by us. The injected electrons is assumed to be a power-law Q(γ e , t ) = Q 0 (t )(γ e /γ m ) −p for γ e > γ m , where Q 0 is related to the injection rate by N inj = γ max γ m Q(γ e , t )dγ e 10 . Then, there are eight free parameters in total in our calculations, i.e., Γ, γ m , B 0 , p, N inj , q, R 0 and η e (or η p ). Particularly, q = 1 is commonly adopted for all calculations unless explicitly stated since the toroidal magnetic field in the ejecta decreases as R −1 . As we will see, this treatment does not markedly impact our main conclusions. Also, p = 2.8 is commonly adopted since the evolution of the low-energy electron distribution is nearly unaffected by p in fast cooling cases. Therefore, there are still six free parameters left. Below, we choose reasonable values for these parameters and perform a sequence of calculations to represent various physical conditions. With the electron spectra being derived numerically, we can then calculate the corresponding synchrotron radiation spectra according to Appendix B.

Testing Calculations
In this section, we first check the significance of SSC cooling in the evolution of the electron distribution in the 10 γ max is the maximum Lorentz factor of electrons and is given by the approximation γ max 10 8 B 1 G −0.5 (Dai & Lu 1999;Huang et al. 2000). So authors may notice that γ max is evolving with time in the results of some calculations.
testing calculations. Four calculations are performed, which are named in form of "Mi" (i = 1, ..., 4). In M1, we set q = 0, ignore SSC and adiabatic cooling so that this should give a "standard" evolution pattern (dN e /dγ e ∝ γ −2 e ) for electrons under synchrotron cooling only. The values for other parameters (shown in Table 1) are taken as the same as those in Uhm & Zhang (2014) in order to compare the results directly. In M2, we set q = 1 to achieve the similar results with the decaying B case shown in Uhm & Zhang (2014). In M3, we introduce the SSC cooling process and compare the results with M2. This should be more close to the realistic situation. At last, we set q = 0 to ignore the effect of decaying B in M4, where the result can clearly show the role played by SSC cooling. The resultant electron distributions and radiation spectra are shown in Figure  1 and Figure 2 respectively for these four calculations. Moreover, we present the cooling rates of different radiation mechanisms in Figure 3. It is not surprising that the indices of the low-energy electron spectra are always strictly −2 for γ e < γ m in M1, and the indices in M2 turn harder along with B 's decreasing as proposed by Uhm & Zhang (2014). In M3 and M4, it is interesting to find that the indices of the low-energy electron spectra are approaching −1 with the elapsing time, and the electrons with γ e < γ m are being cooled mainly via SSC process as shown in Figure 3. The asymptotic value of −1 is consistent with what is predicted theoretically as mentioned. Another natural result is that the electrons are cooled much faster after considering the SSC process. For example, in M3, the minimum Lorentz factor of electrons has already reached 20 at 0.03 s in the observer frame, while it takes 1.5 s for electrons to cool to γ e = 100 in M2. With these testing calculations, we see that SSC cooling can play an important role in determining the electron distribution, at least for in cases considered in previous researches.

Cases in Different Scenarios
Next, we perform numerical calculations by taking the parameters in plausible ranges for GRBs, corresponding to different physical scenarios,

The Poynting-flux-dominated Jet
Let us first consider the case that the jet is a Poyntingflux dominated. According to Equation (16), we can have the upper limit of γ m (η e = 1) by taking the following set of typical parameters: F ν obs = 1 mJy, E peak = 500 keV, z = 1 and Y = 0. In this scenario, we only use the upper limit of γ m in the calculations. One may note that η e = 1 means that L e = L B , which seems to be contrary to the fact of a Poynting-flux-dominated jet. However, we would see that the corresponding results are representative. In addition, we perform two groups of calculations, in which R 0 is taken to be 10 15 cm (called Group PJR15) and 10 14 cm (called Group PJR14), respectively. Then, in each group, we assume a series of values for Γ, and the corresponding B 0 can be obtained via Equation (11). For R 0 = 10 15 cm, we perform five calculations named in form of "PJR15ΓN" with N denoting the value of Γ, which ranges from 50 to 10 3 (see details in Table 2). For R 0 = 10 14 cm, we perform four calculations named in form of "PJR14ΓN", in which Γ ranges from 10 2 to 10 3 (see details in Table 3).
From the results of PJR15Γ1000, we notice that the electron cooling is dominated by the synchrotron radiation for γ e > 10 2 , whereas the adiabatic expansion dominates the cooling of electrons with γ e < 10 2 (see Figure 4). It indicates that SSC cooling is not significant in this case, which also holds for other cases with even greater B 0 in Group PJR15. On the other hand, calculating the SSC cooling rates is extremely time-consuming since it is difficult to achieve the integral convergence in our code. So, we do not include the SSC cooling effect in this group of calculations. The corresponding results for electron distributions, flux spectra, and cooling rates are shown in Figures 5, 6 and 7, respectively. In the results of PJR15Γ1000 and PJR15Γ600, it is seen that along with the decaying of B , adiabatic cooling becomes more and more dominant for low-energy electrons sincė γ e,syn ∝ B 2 ∝ R −2 andγ e,adi ∝ R −1 . The decreasing ofγ e,adi along R leads to an electron spectrum that is harder than −1, just like what is done by the decreasinġ γ e,syn (see results of M2). The low-energy indices of the flux density spectra in PJR15Γ1000 (α, F ν ∝ ν α ) can be larger than 0 and α in PJR15Γ600 are within [-0.5,0]. In results of PJR15Γ300, PJR15Γ100 and PJR15Γ50, the cooling of all electrons is dominated by the synchrotron radiation only. In Figure 5, the decaying of B leads the electron spectra to become harder than -2 as expected in PJR15Γ300. However, this does not occur in PJR15Γ100 or PJR15Γ50. Comparing with the results of Uhm & Zhang (2014), the electron spectrum seems to become harder when the electrons are cooled at a decreasing cooling rate. A history/experience of decreasing cooling rate for an electron is the key factor to result in a hard spectrum. For relatively large B and R, the cooling timescale of electrons is significantly smaller than the dynamical timescale. As a consequence, the history of the decreasing cooling rate of an electron is too short to take effect.
In calculations of Group PJR14, the results for electron distributions, flux spectra, and cooling rates are shown in Figures 8, 9 and 10, respectively. In the results of PJR14Γ1000, it can be seen that adiabatic cooling is dominant for low-energy electrons due to the relatively large Γ and small B , R 0 . The indices of the low-energy electron spectra in PJR14Γ1000 are harder than −1 and α also becomes harder than 0 at frequency ∼ 1 keV when t obs > 0.5 s. In PJR14Γ600, synchrotron cooling is dominant at early times, while adiabatic cooling becomes dominant at late stages. The low-energy electron spectra and the flux spectra are harder than the standard ones, but not so hard as those in PJR14Γ1000. In the results of PJR14Γ300 and PJR14Γ100, we can find that synchrotron cooling is always dominant and the low-energy electron spectra and the flux spectra are just similar to the standard ones. Here again, in PJR14Γ300 and PJR14Γ100, the cooling timescale of an electron under relatively large B is significantly smaller than the dynamical timescale, and the mechanism of "decreasing synchrotron cooling rate" cannot work to make the electron spectra significantly harder than −2.
For Poynting-flux-dominated jets considered here, we confirm that the decreasing synchrotron cooling rate (decreasing B ) will lead to hard electron spectra and flux spectra, according to the results of our calculation Groups PJR15 and PJR14. The γ m used in PJR15 and PJR14 are one or two orders smaller than those used in (Uhm & Zhang 2014), to meet the physical condition of L B ≥ L e . Although γ m used here is only the upper limit, the results from Calculations PJR15 and PJR14 should be representative for three kinds of cooling patterns in this scenario, i.e., hard electron spectra caused by decreasingγ e,adi , or decreasingγ e,syn , and normal spectra under large B . Moreover, we notice that B should not be too large in this scenario, otherwise the synchrotron cooling timescale would be much shorter than the dynamical timescale and the effect of decaying B is weakened. This is also the reason that we did not explore the regime of L B L e (when SSC cooling is not important). In other words, if we want to use the decreasingγ e,syn to work for hard electron spectra, there should be a lower limit ζ for the extent of fast cooling, i.e., ζ ≤ t c /t d < 1. The "unsuccessful" results from PJR15Γ100-PJR15Γ50 and PJR14Γ300-PJR14Γ100 together indicate that ζ should be larger than 10 −5 . Since , the proper range of B for decayingγ e,syn to work is likely to be within roughly 2 orders of magnitude. If this mechanism is true for GRB spectral hardening, the narrow range of B indicates a potential way to probe the magnetic strength of the GRB jet using its spectral characteristics.

The Internal Shock Scenario
Now, we turn to the internal shock model. In this scenario, since we still have little knowledge about the emission radii of GRBs and the crucial parameter η p , we perform four groups of calculations to try to cover various possibilities. For R 0 = 10 15 cm, we consider two situations, i.e., η p = 20 11 or η p = 100. For η p = 20, we show the allowed parameter region in the γ m − Γ diagram (see Figure 11) by combing the restrictions given in Equations (13) and (17). In this group of calculations, called Group IS20R15, four calculations are performed and the positions of the corresponding parameters in the parameter space are marked as black stars in Figure 11. Detailed parameter values in each calculation are listed in Table  4 and each calculation is named in form of "IS20R15N", with N denoting the value of Γ. Similarly, in the calculations of Group IS100R15 (η p = 100), we perform thirteen calculations of which the corresponding information can be seen in Table 5 and Figure 12. These thirteen calculations are classified into four subgroups according to different values used for γ m . Each calculation is named in form of "IS20R15WΓN" with letter "W" (e.g., A, B etc.) distinguishing different subgroups and N denoting the value of Γ. When R 0 = 10 14 cm is adopted, we also consider two situations, i.e., η p = 10 or η p = 100. Six calculations are performed in Group IS10R14 (see Table  6 and Figure 13), while ten calculations are performed in Group IS100R14 (see Table 7 and Figure 14).
In the calculations of Group IS20R15, we notice that γ m = 1.5 × 10 4 and L e ≈ L B . The results for electron distributions, flux spectra, and cooling rates are shown in Figures 15, 16, and 17 respectively. From these results, we find that the electron spectra and the flux spectra can be hard enough to match the observations only when adiabatic cooling is dominant (in IS20R15Γ1300). For other three cases (IS20R15Γ430-IS20R15Γ86),γ e,SSC is only comparable toγ e,syn for electrons of γ e ≤ 10 2 , and the indices of the low-energy electron spectra are slightly harder than −2.
In the calculations of Group IS100R15, we can see γ m ranges from 7 × 10 3 to 10 5 and L e /L B ranges from 0.08 to 3000. According to the results (see Figures  18, 19 and 20), we find: 1, for Subgroup IS100R15A (L e /L B = 3000),γ e,SSC is always dominant for electrons of γ e < γ m and the resulting indices of electron spectra are ∼ −1.3; 2, for Subgroup IS100R15B (L e /L B = 200),γ e,SSC is dominant for electrons of γ e < 10 4 and the resulting indices of electron spectra are ∼ −1.4, except for IS100R15BΓ460, in whichγ e,adi becomes dominant at the late time and the low-energy electron indices can be even harder than −1; 3, for Subgroup IS100R15C (L e /L B = 0.3),γ e,syn is always larger thaṅ γ e,SSC . The electron spectra resemble the standard ones in IS100R15CΓ580 and IS100R15CΓ58, while the electron indices in IS100R15CΓ1900 are becoming harder than −1 due to the dominance of adiabatic cooling; 4, for Subgroup IS100R15D (L e /L B = 0.08),γ e,syn is always larger thanγ e,SSC . The electron spectra resemble the standard ones in IS100R15DΓ1200 and IS100R15DΓ120, while the electron indices in IS100R15DΓ3900 are becoming even harder than 0 due to the dominance of adiabatic cooling.
In the calculations of Group IS10R14, we find (see Figures 21, 22 and 23): 1, for Subgroup IS10R14A (L e /L B = 30),γ e,SSC is larger thanγ e,syn for electrons of γ e < 3 × 10 3 , and the resulting indices of electron spectra are ∼ −1.5 for IS10R14AΓ580 and IS10R14AΓ58. For IS10R14AΓ1900,γ e,adi is always dominant for electrons of γ e < γ m and the electron indices are becoming harder than 0; 2, for Subgroup IS10R14B (L e /L B = 2),γ e,SSC is slightly larger thanγ e,syn for electrons of γ e < 3 × 10 2 , and the resulting indices of electron spectra are ∼ −1.8 for IS10R14BΓ770 and IS10R14BΓ230. For IS10R14BΓ2300,γ e,adi is becoming increasingly dominant for γ e < γ m and the electron indices are becoming harder than 0.
In the calculations of Group IS100R14, γ m ranges from 5 × 10 3 to 10 5 and L e /L B ranges from 2 to 3 × 10 5 . According to the results (see Figures 24, 25 and 26), we find: 1, for Subgroups IS100R14A and IS100R14B (L e /L B is 3 × 10 5 or 8.5 × 10 3 ),γ e,SSC is always much larger than bothγ e,syn andγ e,adi for electrons of γ e < γ m , the resulting indices of electron spectra are ∼ −1; 2, for Subgroup IS100R14C (L e /L B = 30),γ e,SSC is larger than bothγ e,syn andγ e,adi for γ e < 3 × 10 3 , the resulting indices of electron spectra are ∼ −1.5. 3, for Subgroup IS100R14D (L e /L B = 2), sinceγ e,adi is becoming increasingly dominant for low-energy electrons, the indices of electron spectra are getting harder than 0. The results of IS100R14DΓ230 are similar to those of IS10R14BΓ230 since their parameters are actually the same. So we have not shown them as subfigures in the relevant figures of this calculation group.
To sum up, in cases of a large Γ (Γ > 10 3 ), adiabatic cooling is the most dominant process for low-energy electrons' cooling, making the electron indices harder than −1, or even 0. For cases of L e /L B > 1, and when Γ is not too large, SSC cooling dominates over synchrotron cooling, resulting in electron spectra with indices ranging from ∼ −1 to −2 with the decreasing of L e /L B . When L e /L B ≤ 1 is met while Γ is still not too large, synchrotron cooling will then take over and the resulting spectra slightly deviate from the standard ones. These three kinds of cooling patterns are generally consistent with previous numerical results in Bošnjak et al. (2009) and Daigne et al. (2011) 12 . The conditions for the dominance of adiabatic cooling or SSC cooling in these results are consistent with the analyses in Section 4. In general, the combination of a large γ m and small B , R will favor the dominance of adiabatic cooling or SSC cooling, rather than synchrotron cooling according to Equations (21) and (27).
In reality, the observed minimum variability timescales of GRBs can be used to derive the internal shock radii, which typically gives R 0 = 10 14 cm or smaller. We have only explored the region of R 0 ≥ 10 14 cm in our calculations. However, according to Equations (26) and (27), one will realize that a smaller R 0 will enhance the SSC cooling rate due to the increase of the radiation energy density, i.e., U γ ∝ R −2 . Therefore, for a smaller R 0 , SSC cooling will be more significant. Also, we can expect that the results of a smaller R 0 can be well differentiated by values of L e /L B and Γ.

DISCUSSION AND CONCLUSIONS
We have developed a code to solve the continuity equation of electrons in GRB ejecta, and analyzed the roles played by three cooling mechanisms (synchrotron, SSC, adiabatic cooling) and the effect of decaying magnetic field in determining the electron/flux spectra in both the fireball and the Poynting-flux-dominated regimes. By exploring the parameter space and calculating the corresponding electron spectra, flux spectra, and electron cooling rates, we find that the hardening of the electron spectra can be attributed to the synchrotron radiation with a decaying B , or the dominance of adiabatic cooling or the dominance of SSC cooling. Therefore, it is essential to coherently consider them together in future studies of the GRB prompt emission. The numerical method, as proposed in this paper, has the advantage over the analytical method in solving equations involving several factors simultaneously. According to our results, the low-energy spectra of GRBs could be explained by the synchrotron radiation from either a Poynting-fluxdominated jet or an internal shock, although some shortcomings may exist for the two scenarios.
SSC cooling of electrons in the KN regime has been proposed to solve the fast cooling problem previously. In this paper, our analyses confirm that SSC cooling is crucial in the internal shock scenario when we take 12 In our results, spectral indices α are strictly smaller than 0 for cases when SSC cooling dominates. This is consistent with the limit of α > −0.1 given in Barniol Duran et al. (2012)  parameters deduced from typical observational characteristics. A sequence of numerical calculations further reveal that SSC cooling will result in electron spectra with low-energy indices ranging from −2 to ∼ −1. The physical condition for SSC cooling to be dominant is L e /L B > 1. In order to match the observations, i.e., dN e /γ e ∼ γ −1 e , F ν ∼ ν 0 , the condition should be more strict, i.e., L e /L B ≥ 10 4 according to the results from the calculation Groups IS100R15 and IS100R14. This condition means ξ e /ξ B ≥ 10 4 within the internal shock model. It is still hard to understand why the energy fraction in the magnetic field could be so small. However, a possible explanation is that the B that cools the electrons here is much smaller than the magnetic field in the acceleration region. As proposed by Rossi & Rees (2003), Pe'er & Zhang (2006), Lemoine (2013), and Zhao et al. (2014), the magnetic field may accumulate in a small region (but carries the majority of the magnetic energy created by the shock) just near the shock front, and the magnetic field behind the shock decays rapidly with the distance from the front. This possibility is hinted in particle-incell simulations of shocks (e.g., Medvedev et al. 2005;Chang et al. 2008). Since the electrons may be cooled in the downstream region behind the shock, then the B they stream through should be smaller than those near the shock. Moreover, to achieve the electron spectra with low-energy indices of ∼ −1, η p ≥ 10 2 is preferred to take γ m ≥ 10 4 . It indicates that only ∼ 1% electrons are accelerated and they carry ∼ 10% of the energy when they cross the shock. This should correspond to the cases in which the relativistic electron-ion shock is of low magnetizations (σ ≤ 10 −3 ) according to the simulations of Sironi & Spitkovsky (2011) and Sironi et al. (2015).
Our numerical results also support the idea that the hardening of the electron spectra purely by SSC cooling in KN regime can be only up to, but not equal to −1 (Barniol Duran et al. 2012). So the fast cooling problem could not be fully solved by SSC cooling, since we have some GRB spectra with α > 0. However, for this small fraction of GRBs, they may be well explained when adiabatic cooling is the dominant cooling process for lowenergy electrons. This requires that the ejecta is moving at an extremely relativistic speed, i.e., Γ > 10 3 . In the framework of the internal shock model, this condition may be naturally fulfilled since the initial energy of the ejecta is mainly turned into kinetic energy. On the other hand, the lower limit for Γ in three GRBs set by the Fermi team is ∼ 1000 (Abdo et al. 2009a,b,c). Thus it is possible that adiabatic cooling may be dominant in a few GRBs.
For a Poynting-flux-dominated jet (L B ≥ L e ), we confirm that the low-energy indices of the electron spectra can be harder than ∼ −1 when B is under decaying with R (Uhm & Zhang 2014), as shown in the results of PJR15Γ600 and PJR14Γ600. By using the constraints from observations and performing a sequence of calculations (see Section 5.2.1), we further reveal that B should not be too large in this scenario, otherwise the synchrotron cooling timescale will be much shorter than the dynamical timescale and the effect of the decaying B is weakened. This feature provides a way to identify the effect of SSC cooling from the effect of decaying B .
Although spectra hardening could be achieved by two mechanisms, i.e., SSC cooling in KN regime or the effect of decaying B , they may be distinguished in observations. For SSC cooling in KN regime, the hardening of electron spectra is a result of the "intrinsic" scale relation ofγ e,tot γ e,SSC ∝ γ −1 e (Derishev et al. 2001;Bošnjak et al. 2009;Wang et al. 2009;Nakar et al. 2009;Fan 2010;Daigne et al. 2011). It can be seen that the hardening of electron spectra could be established within a rather short timescale, e.g., t obs < 0.1 s via SSC cooling as in the results of Subgroups IS100R14A and IS100R14B. In contrast, decaying B is an "external" way, and hardening of the electron spectra due to this effect needs a longer time, e.g., t obs ≥ 0.5 s according to the results of PJR15Γ600 and PJR14Γ600. Therefore, if a flux spectrum of form F ν ∼ ν 0 is observed in a single, distinguishable, short pulse (t obs < 0.1 s) for a particular GRB, then it is more likely that the SSC mechanism should be the dominant cooling process. On the other hand, spectral lags and E peak evolution patterns are found to be related with broad pulses (with durations of seconds, or called slow component, Gao et al. 2012) rather than quick variabilities (Zhang & Yan 2011;Uhm & Zhang 2016a). According to R ∼ Γ 2 cδt obs , the emission radius is large for large δt obs , which would weaken the effect of SSC cooling. The large emission radius is also contrived within the internal shock model. Thus, for these broad pulses, the effect of decaying B should be more important if the evolution of spectral hardening (and softening) could be well matched with numerical results.
No general consensus on GRB jet properties (e.g., jet composition, emission radius) has been reached in the community. The information from the time-resolved lowenergy spectra can help to probe the physical characteristics of the GRB ejecta via our numerical results. As mentioned above, SSC cooling in KN regime works in the scenario of internal shocks (baryon-dominated jet), while the effect of decaying B mainly happens in the scenario of Poynting-flux-dominated jet. Once the timeresolved low-energy spectra hardening is affirmed to be due to a specific mechanism, the jet composition could also be inferred simultaneously. Furthermore, an overall comparison of the results among the four calculation groups within the internal shock model indicates that α ∼ 0 is more likely to be achieved at a small emission radius (≤ 10 14 cm). This value is well consistent with the results from other methods (Rees & Meszaros 1994;Gupta & Zhang 2008;Kumar & Zhang 2015).
There is another way to solve the fast synchrotron cooling problem. Except for the cooling processes of electrons, their acceleration processes should also be important in determining the final electron distribution. A hard electron energy distribution with index ∼ −1 may be produced by a slow heating process (e.g., Ghisellini & Celotti 1999;Stern & Poutanen 2004;Asano & Terasawa 2009). This possibility was confirmed by Xu & Zhang (2017) recently, who considered the second-order Fermi acceleration in the turbulent reconnection (Zhang & Yan 2011). In future studies, we will incorporate the acceleration term into Equation (5) and investigate the effect in more details.
Some parameters are fixed in our calculations, such as N inj , Γ etc. Although these parameters may be variable in reality, our results can still be a useful baseline for further detailed explorations. Actually, the variation of these parameters can be easily taken into account in our numerical code, by setting the boundary conditions as time-dependent functions. Furthermore, more physical processes could be considered instantly. For example, the scattering of electrons by various external radiation fields (Yan et al. 2016) may also be an important cooling mechanism. Additionally, the electrons that are not efficiently accelerated (the low-energy electrons of Maxwellian distribution) may provide considerable lowenergy seed photons. For the complete test of a GRB model, it is necessary to consider the dynamics of the GRB jet at much earlier stages. For example, a Poyntingflux-dominated jet may undergo accelerating when it is emitting (Uhm & Zhang 2016b), while the final Lorentz factor of merged ejecta is determined by the momentum of early shells. In the future, we will improve our code to simulate the dynamical processes of GRBs more physically, deriving the light curves and spectra at the same time. Such complete theoretical fitting to the observed spectra and light curves will give more clues on the characteristics of the emission region. We note that the Hard X-ray Modulation Telescope (HXMT) (Li 2007;Xie et al. 2015) launched by China recently can cover an energy range of 1 -250 keV. It will be efficient in collecting the time-resolved low-energy spectra of GRBs in the near future. Our methodology proposed here will be helpful to make use of the HXMT data.
We thank the anonymous referee for valuable suggestions. We also thank Z. Lucas Uhm, Mi-Xiang Lan, Liang Li and Bing Li for helpful discussion. J.J.G thanks Pengfei Chen for useful help on the numerical methods. This work is partially supported by the National Postdoctoral Program for Innovative Talents

NUMERICAL METHOD
In this appendix, we present the discretization procedure to solve Equation (5), i.e., the CIP method (see Yabe et al. 2001 for a general review). For an one-dimensional non-linear equation it is expedient to separate the solution procedure into two fractional steps. One is the advection phase, and the other is the non-convection phase, wheref = ∂f /∂x,Ĝ = ∂G/∂x stands for the spatial derivative of f and G.f should be solved together with f in the CIP method, which is crucial to obtain the propagation of the spatial derivative during the evolution. If we assume that the profile between two adjacent points can be interpolated by the cubic polynomial F (x) = ax 3 + bx 2 + cx + d, then the solution at grid i can be evolved from step n to step n + 1 by where , ξ = −u∆t, and ∆t is the time step. Note that if u < 0 (just the case for electron cooling), one should replace i − 1 with i + 1, and ∆x with −∆x. f * i andf * i in Equations (A6) and (A7) are the the intermediate solution from the non-convection phase, and can be obtained through centered finite-difference of (A4) and (A5), i.e., For the specific case in this paper, since the thermal Lorentz factor γ e ranges from 10 to 10 7 , it is necessary to solve Equation (5) in the logarithm space of γ e . We determine log 10 γ e = x, so that dNe dx = ln 10γ e dNe dγ e , and the Equation (5) where dNe dx is actually what we want to solve in our code. In all the calculations, we set the range of x to be [1,8], and the total grids number N num to be 401. The time step for every evolution is determined by the Courant condition where γ max is the maximum thermal Lorentz factor and can be given by the approximation γ max 10 8 B 1 G −0.5 (Dai & Lu 1999;Huang et al. 2000). However, the uncertainty of the real value of γ max has little influence on the evolution of the electron distribution and the spectra we focus on.

FORMULATIONS FOR RADIATION
In the co-moving frame, the synchrotron radiation power at frequency ν is (Rybicki & Lightman 1979) where ν c = 3q e B γ 2 e /(4πm e c), F (x) = x +∞ x K 5/3 (k)dk, and K 5/3 (k) is the Bessel function. The synchrotron seed photon spectra can then be calculated as (Fan et al. 2008) where T ≈ ∆/c is the time that the synchrotron radiation photons stay within the ejecta, n (γ e ) = dNe/dγ e 4π∆R 2 is the co-moving electron number density, and ∆ ≈ R/Γ is the co-moving width of the ejecta. One should notice that ∆ does not appear in our calculations since the ∆ in T and n (γ e ) are canceled out. Here, for simplicity, we have considered only the single scattering case for SSC cooling and ignored the multiple scattering process.
If we ignore the effect of the equal-arrival-time surface (EATS, Waxman 1997;Granot et al. 1999;Huang et al. 2007;Geng et al. 2017), then the observed spectral flux can be expressed as where ν = (1 + z)ν obs /D, and D = 1/[Γ(1 − β cos θ)] is the Doppler factor. In this work, the luminosity distance D L is obtained by adopting a flat ΛCDM universe, in which H 0 = 71 km s −1 , Ω m = 0.27, and Ω Λ = 0.73. For all the calculations in this work, the burst is assumed to be at a cosmological redshift z = 1. If we take the EATS effect into account, the observed spectral flux should be where θ j is the half-opening angle of the jet. The integration of θ is performed over an elliptical surface (or a sequence of R θ ), which is determined by (Geng et al. 2016) from which R θ can be derived for a given θ. However, we found that there is little difference for the spectrum calculated from the case with the EATS effect and the case without in our calculations. Also the EATS effect has no significant influence on the low-energy indices of the flux spectra. So we just show the spectra calculated without the EATS effect in this paper. Note. -In this group of calculations, one should note that R0 = 10 15 cm is adopted, whereas the starting radius at which the jet begin to produce emission is Rs = 10 14 cm. This is only to achieve the same initial conditions as those in Uhm & Zhang (2014). However, Rs = R0 is commonly used all through this paper.         Table 1). The orange dashed lines are the standard fast cooling pattern, i.e., dNe/dγ e ∝ γ −2 e , and the grey dashed lines present the expected cooling pattern according to the observations, i.e., dNe/dγ e ∝ γ −1 e . The epochs shown for each case are different since the electron cooling timescales are different. The lower panel of each case is the negative spectral index of the electron spectrum, i.e., χ(γ e ) = dNe/dγ e . In the standard calculation, M1, the electron spectrum shows a typical broken power-law profile, with the spectral indices being −(p + 1) above γ m and being -2 below γ m . In other cases, the spectral indices below γ m are notably harder than -2. In the standard calculation, M1, it gives the typical fast cooling spectrum Fν ∝ ν −1/2 below νm (νm ∝ γ 2 m ). For even lower frequency, the spectrum is Fν ∝ ν 1/3 , which is the profile of the low frequency part of a single electron's synchrotron spectrum. In other cases, the spectral indices below νm can be harder than −1/2 and can approach 0.    Table 2).   Table 3).
The parameter space for the internal shock model when R 0 = 10 15 cm and ηp = 20 are adopted. The shadow region is the plausible region for the γ m -Γ couple constrained from Equations (13) and (19). The magenta dashed line presents the condition for fast cooling, while the grey dashed and the red dashed lines show the lower limit (γ m < ηpmp/me) and the upper limit (deduced by ξ B ≤ 1 with Equation (17)) of γ m , respectively. Γ is assumed to range from 10 to 10 4 according to previous researches on GRB jets (Liang et al. 2010(Liang et al. , 2013. The solid lines are the relationship between γ m and Γ using Equation (11), when different values for B (denoted by different colours) are adopted and E peak is set to be 500 keV typically. Particularly, the cyan dotted line represents the minimum value of B for γ m -Γ couples to overlap with the shadow region. The positions of the parameters in the four calculations of Group IS20R15 (see Table 4) are marked by star symbols.  Fig. 12.-The parameter space for the internal shock model when R 0 = 10 15 cm and ηp = 100 are adopted. The meanings of lines are similar to those explained in Figure 11. The positions of the parameters in the thirteen calculations of Group IS100R15 (see Table 5) are marked by star symbols.
The parameter space for the internal shock model when R 0 = 10 14 cm and ηp = 10 are adopted. The meanings of lines are similar to those explained in Figure 11. The positions of the parameters in the six calculations of Group IS10R14 (see Table 6) are marked by star symbols.
14.-The parameter space for the internal shock model when R 0 = 10 14 cm and ηp = 100 are adopted. The meanings of lines are similar to those explained in Figure 11. The positions of the parameters in the ten calculations of Group IS100R14 (see Table 7) are marked by star symbols.  Table 4).   Table 5).       Table 6).    -The evolution of the electron energy spectrum for the ten cases in Group IS100R14 (see Table 7).