© The Authors 2026

1. Introduction

Gamma-ray bursts (GRBs) are among the most energetic explosions in the universe, and understanding their central engine mechanisms is crucial to probing jet dynamics and radiation processes. Currently, the central engines of GRBs are commonly attributed to two models: a magnetar (Duncan & Thompson 1992; Dai & Lu 1998a,b; Metzger 2010; Liu et al. 2015) and a hyperaccreting stellar-mass black hole (Blandford & Znajek 1977; Paczynski 1991; Woosley 1993; Popham et al. 1999; Narayan et al. 2001; Lei et al. 2013; Xie et al. 2016; Lei et al. 2017). In the magnetar scenario, a rapidly rotating, strongly magnetised neutron star powers relativistic jets through accretion; continued energy injection may produce an extended X-ray plateau – a feature typically observed in bursts with lower energetics and longer plateau durations. Meanwhile, hyperaccreting black holes originate from either massive stellar collapse or compact binary mergers. These systems generate relativistic jets via the Blandford–Znajek (BZ) mechanism (Blandford & Znajek 1977; Lee et al. 2000a,b; McKinney 2005; Lei et al. 2013; Liu et al. 2015; Lei et al. 2017), driven by magnetic fields or through the annihilation of neutrino–antineutrino pairs (NDAF) (Ruffert et al. 1997; Popham et al. 1999; Di Matteo et al. 2002; Liu et al. 2007; Zalamea & Beloborodov 2011; Xue et al. 2013; Song et al. 2016; Xie et al. 2016).

Some GRBs are attributed to the magnetar model, whereas the hyperaccreting black hole model can explain the majority of GRBs, particularly high-luminosity and ultra-long GRBs. Here, we focus on the hyperaccreting black hole system. The true jet luminosity, inferred by excluding X-ray emission from the central engine, can be linked to the NDAF and BZ mechanisms using parameters such as the isotropic radiated energy E_γ,iso, the isotropic kinetic energy E_k,iso, redshift z, and the jet half-opening angles θ_j. Fan & Wei (2011) estimated the accretion disc mass for short GRBs to range from 0.01 M_⊙ to 0.1 M_⊙, suggesting that the binary neutron star merger scenario can plausibly account for their observed characteristics. Liu et al. (2015) and Song et al. (2015, 2016) further pointed out that the NDAF mechanism faces difficulties in explaining high-luminosity and ultra-long GRBs, particularly those requiring large accretion disc masses. Yi et al. (2017) found that the BZ mechanism matches the selected observational data better than the NDAF mechanism, a result further supported by Xue et al. (2013). Du et al. (2021) refined the jet-accretion rate constraints for high-luminosity GRBs within the BZ scenario.

Although previous studies have primarily focused on energy estimates or statistical analyses, some investigations have explored how light-curve variability reflects the central engine mechanism. The morphology of the light curve provides crucial insights into this distinction. Lei et al. (2007) linked complex light-curve structures to BZ-driven jet precession and nutation, and successfully fitted the light curves of five GRBs. Kumar et al. (2008a) explained the prompt emission phase and afterglow plateau through fallback accretion’s impact on jet luminosity. Lei et al. (2017) provided a general luminosity fitting formula for both NDAF and BZ mechanisms, analysing light-curve evolution across prompt and afterglow phases. Thus, distinguishing the central engine based on the instantaneous light-curve features holds significant research value.

In this work, we simulated light curves during the prompt radiation phase for different initial parameters of hyperaccreting black hole systems to compare the light curves for BZ and NDAF scenarios. To validate our theoretical models, we selected GRBs exhibiting single-peaked FRED profiles observed by Swift and compared their theoretical light curves. Through this comparison, we establish a new link between central engine physics and light-curve morphology. Section 2 describes the hyperaccreting black hole and prompt emission mechanisms. Section 3 details the sample selection and processing methods. Section 4 presents the results, with Sections 5 and 6 discussing and concluding the study.

2. Hyperaccreting black hole

Stellar-mass black hole hyperaccretion systems are widely considered as prominent candidate models for the central engine of GRBs (Narayan et al. 1992; MacFadyen & Woosley 1999). These systems primarily involve two candidate mechanisms: the magnetically driven BZ mechanism and the NDAF mechanism. This section reviews both mechanisms and compares their effects on the light curves of GRBs.

2.1. BZ mechanism

The BZ mechanism, first proposed by Blandford & Znajek (1977), was initially introduced to explain the energy generation in active galactic nuclei. Later, Lee et al. (2000a) and Lee et al. (2000b) extended this mechanism to explain the central engine of GRBs. They suggested that when a rotating Kerr black hole is threaded by a large-scale multipolar magnetic field from the accretion disc, the magnetic field efficiently extracts the black hole’s rotational energy, driving relativistic jets. In this process, the magnetic field lines, dragged by the black hole’s spin, induce an electric field and current near the event horizon, forming a stable magnetosphere. The current generated exerts an Ampère force on the black hole, causing it to decelerate. Consequently, the black hole’s rotational energy is transmitted outward as Poynting flux, which drives the relativistic jets.

In the BZ mechanism, the energy extraction rate from the rotating black hole’s magnetic field can be approximated as (Lei et al. 2017)

$$\begin{array}{c}\hfill L_{\mathrm{BZ}}=1.7\times {10}^{20}{a}_{\bullet }^2{m}_{\bullet }^2{B}_{\bullet ,G}^{2}F(a_{\bullet })\text{erg}{\text{s}}^{-1},\end{array}$$(1)

where m_∙ = M_∙/M_⊙ is the dimensionless black hole mass; B_∙ refers to the magnetic field strength at the black hole horizon, B_∙,G = B_∙/1 G; and F(a_∙) is a dimensionless function of the black hole’s spin parameter a_∙, given by $F(a_{\bullet })=[(1+q^2)/q][(q+1/q)arctan(q)-1]$, where q is defined as $q=a_{\bullet }/(1+\sqrt{1-a_{\bullet }^2})$.

By assuming equilibrium between the accretion disc’s gas pressure and the magnetic pressure at the disc’s inner radius (Mészáros & Rees 1997), we obtain

$$\begin{array}{c}\hfill \frac{B_{\bullet }^2}{8\pi }=P_{\text{in}}\sim \frac{\dot{M}c}{4\pi r_{\mathrm{H}}^2},\end{array}$$(2)

where Ṁ represents the accretion rate at the inner edge of the accretion disc, and r_H is the radius of the black hole’s event horizon, given by $r_{\mathrm{H}}=(1+\sqrt{1-a_{\bullet }^2})r_{\mathrm{g}}$. We then obtain

$$\begin{array}{c}\hfill B_{\bullet }\simeq 7.4\times {10}^{16}{\dot{m}}^{1/2}{m}_{\bullet }^{-1}{(1+\sqrt{1-a_{\bullet }^2})}^{-1}\text{G}.\end{array}$$(3)

Finally, the luminosity of the BZ mechanism can be approximated as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill L_{\mathrm{BZ}}& =9\times {10}^{53}{a}_{\bullet }^2\dot{m}X(a_{\ast })\text{erg}{\text{s}}^{-1}\hfill \\ \hfill & \simeq 1.5\times {10}^{53}{a}_{\bullet }^2\dot{m}\text{erg}{\text{s}}^{-1},\hfill \end{array}\end{array}$$(4)

where ṁ = Ṁ/(M_⊙ s⁻¹ is the dimensionless mass-accretion rate of a black hole.

To establish a connection between the BZ power and the observed gamma-ray luminosity, we use the following relation:

$$\begin{array}{c}\hfill \eta L_{\mathrm{BZ}}=f_{\mathrm{b}}{L}_{\gamma ,\text{iso}}=L_{\gamma },\end{array}$$(5)

where f_b is the beaming factor of the jet, and η is the efficiency of converting BZ power into gamma-ray radiation. For this study, we adopt typical values of f_b = 0.01 and η = 0.2.

In the BZ mechanism, considering the energy and angular momentum conservation (Lee et al. 2000a,b), the evolution of the black hole mass and angular momentum can be expressed as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{M}_{\bullet }{c}^2}{\mathrm{d}t}={\dot{M}}_{\bullet }{c}^2{E}_{\text{ms}}-L_{\mathrm{BZ}},\end{array}$$(6)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{J}_{\bullet }}{\mathrm{d}t}={\dot{M}}_{\bullet }{L}_{\text{ms}}-T_{\mathrm{B}},\end{array}$$(7)

where T_B represents the torque produced by the interaction of the black hole with the accretion disc, and can be written as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill T_{\mathrm{B}}& =\frac{L_{\mathrm{BZ}}}{{\mathrm{\Omega }}_{\mathrm{F}}}=1.8\times {10}^{49}{a}_{\bullet }\dot{m}{m}_{\bullet }F(a_{\ast })\text{g}{\text{cm}}^2{\text{s}}^{-2}\hfill \\ \hfill & \simeq 1.2\times {10}^{49}{a}_{\bullet }\dot{m}{m}_{\bullet }\text{g}{\text{cm}}^2{\text{s}}^{-2},\hfill \end{array}\end{array}$$(8)

where Ω_F = 0.5 Ω_∙ (Lee & Kim 2000; Lee et al. 2000b), and ${\mathrm{\Omega }}_{\bullet }=\frac{a_{\bullet }c}{2r}=\frac{a_{\bullet }{c}^3}{2G{M}_{\bullet }(1+\sqrt{1-a_{\bullet }^2})}$.

Here, E_ms and L_ms correspond to the specific energy and specific angular momentum at the maximum stable circular orbit radius r_ms (Novikov & Thorne 1973), and are given by

$$\begin{array}{c}\hfill E_{\text{ms}}=\frac{4\sqrt{R_{\text{ms}}}-3a_{\bullet }}{\sqrt{3R_{\text{ms}}}},\end{array}$$(9)

$$\begin{array}{c}\hfill L_{\text{ms}}=\frac{G{M}_{\bullet }}{c}\frac{2(3\sqrt{R_{\text{ms}}}-2a_{\bullet })}{\sqrt{3R_{\text{ms}}}},\end{array}$$(10)

where R_ms = r_ms/r_g, and the expression for r_ms is given by

$$\begin{array}{c}\hfill r_{\text{ms}}=r_{\mathrm{g}}[3+Z_2-\text{sgn}(a_{\bullet }){\left([3-Z_1](3+Z_1+2Z_2)\right)}^{1/2}].\end{array}$$(11)

For 0 ≤ a_∙ ≤ 1, the quantity Z₁ is given by $Z_1\equiv 1+{(1-a_{\bullet }^2)}^{1/3}[{(1+a_{\bullet })}^{1/3}+{(1-a_{\bullet })}^{1/3}]$, while Z₂ is defined as $Z_2\equiv {(3a_{\bullet }^2+Z_1^2)}^{1/2}$ (Bardeen et al. 1972; Novikov 1998; Kato et al. 2008).

Building upon this foundation and incorporating the definition of the intrinsic angular momentum, one can derive the expression for the spin as $a_{\bullet }=\frac{J_{\bullet }c}{G{M}_{\bullet }^2}$. The corresponding evolution equation for the Kerr black hole can then be derived as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{a}_{\bullet }}{\mathrm{d}t}=\frac{(\dot{M}{L}_{\text{ms}}-T_{\mathrm{B}})c}{G{M}_{\bullet }^2}-2a_{\bullet }\left(\frac{\dot{M}{c}^2{E}_{\text{ms}}-L_{\mathrm{BZ}}}{M_{\bullet }{c}^2}\right).\end{array}$$(12)

This equation governs the time evolution of the black hole’s spin, derived by considering the BZ mechanism in the context of energy extraction from the black hole.

2.2. NDAF model

The mechanism of neutrino annihilation is an important channel of energy release, first proposed by Popham et al. (1999) to explain GRBs. The typical accretion rate associated with GRBs is in the range of 0.01–10 M_⊙ s⁻¹, resulting in an accretion disc with extremely high temperature and density. Under such conditions, the optical depth increases, photons become trapped, and the viscously dissipated energy can be carried away only via neutrino emission; this constitutes the NDAF mechanism. In this mechanism, a fraction of neutrinos annihilate into electron–positron pairs ($\nu +\overline{\nu }\to e^{-}+e^{+}$), which subsequently annihilate into gamma-ray photons that help launch relativistic jets. Because this process does not require strong magnetic fields, it is regarded as a complementary GRB launching mechanism to the BZ process.

To accurately describe the relativistic accretion system, we adopt a steady-state disc model around a Kerr black hole, incorporating neutrino radiation and rotation within the framework of general relativity, as proposed by Chen & Beloborodov (2007) and Lei et al. (2017). For simplicity, we assume that the NDAF accretion disc is a standard thin disc and that the rotation of the accretion disc follows the Keplerian disc model.

Under the Kerr metric, the basic equations of the NDAF model are as follows:

1. Mass conservation

   $$\begin{array}{c}\hfill \dot{M}=-4\pi r{v}_r\rho H.\end{array}$$(13)

2. Angular momentum conservation

   $$\begin{array}{c}\hfill \dot{M}\sqrt{G{M}_{\bullet }r}\frac{D}{A}=4\pi r^{2}H\alpha P\sqrt{\frac{A}{\mathit{BC}}},\end{array}$$(14)

   where A, B, C, and D are general relativistic correction factors, given as

   $$\begin{array}{c}\hfill A=1-\frac{2G{M}_{\bullet }}{c^{2}r}+{\left(\frac{G{M}_{\bullet }·a_{\bullet }}{c^{2}r}\right)}^2,\end{array}$$

   $$\begin{array}{c}\hfill B=1-\frac{3G{M}_{\bullet }}{c^{2}r}+2a_{\bullet }{\left(\frac{G{M}_{\bullet }}{c^{2}r}\right)}^{3/2},\end{array}$$

   $$\begin{array}{c}\hfill C=1-4a_{\bullet }{\left(\frac{G{M}_{\bullet }}{c^{2}r}\right)}^{3/2}+3{\left(\frac{G{M}_{\bullet }·a_{\bullet }}{c^{2}r}\right)}^2,\end{array}$$

   $$\begin{array}{c}\hfill D=Bf.\end{array}$$

3. Energy conservation

   $$\begin{array}{c}\hfill Q^{+}=Q^{-},\end{array}$$(15)

   where Q⁺ = Q_vis and Q⁻ = Q_ν + Q_photo + Q_adv. Here, Q_vis is the viscous heating rate, and the cooling terms (Q⁻) include neutrino cooling Q_ν, photon emission Q_photo, and advection energy Q_adv.

4. The total pressure on the NDAF disc

   $$\begin{array}{c}\hfill \begin{array}{cc}\hfill P& =\frac{11}{12}a{T}^4+\frac{\rho kT}{m_{\mathrm{p}}}(1+\frac{3X_{\text{nuc}}}{4})+\frac{2\pi hc}{3}{\left(\frac{3}{8\pi m_{\mathrm{p}}}\right)}^{4/3}\hfill \\ \hfill & \times {\left(\frac{\rho }{{\mu }_e}\right)}^{4/3}+\frac{u_{\nu }}{3}.\hfill \end{array}\end{array}$$(16)

The total energy emitted by neutrinos from the NDAF accretion disc per unit time represents the neutrino luminosity. The neutrino radiation luminosity of the NDAF is given by

$$\begin{array}{c}\hfill L_{\nu }=4\pi {\displaystyle {\int }_{r_{\text{ms}}}^{r_{\text{out}}}}{Q}_{\nu }r\mathrm{d}r,\end{array}$$(17)

where r_out represents the outer radius of the accretion disc and Q_ν denotes the neutrino cooling rate. The neutrino radiation luminosity is obtained by integrating Q_ν over the area.

For the neutrino annihilation luminosity, Ruffert et al. (1997), Popham et al. (1999), Rosswog et al. (2003) and others proposed dividing the NDAF disc into many grid points. Neutrinos are denoted as k, while antineutrinos are labelled as k′. By summing and integrating over all the grid points, the neutrino annihilation luminosity at each point is obtained:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\overline{E}}_{\nu \overline{\nu }}& =A_1{\displaystyle \sum _k}\frac{l_{{\nu }_i}^k}{d_k^2}{\displaystyle \sum _{k\prime }}\frac{l_{{\nu }_i}^{k\prime }}{d_{k^{\prime }}^2}({ϵ}_{{\nu }_i}^k+{ϵ}_{{\overline{\nu }}_i}^{k\prime }){(1-cos{\theta }_{kk\prime })}^2\hfill \\ \hfill & +A_2{\displaystyle \sum _k}\frac{l_{{\nu }_i}^k}{d_k^2}{\displaystyle \sum _{k\prime }}\frac{l_{{\nu }_i}^{k\prime }}{d_{k\prime }^2}\left(\frac{{ϵ}_{{\nu }_i}^k+{ϵ}_{{\overline{\nu }}_i}^{k\prime }}{{ϵ}_{{\nu }_i}^k{ϵ}_{{\overline{\nu }}_i}^{k\prime }}(1-cos{\theta }_{kk\prime })\right),\hfill \end{array}\end{array}$$(18)

where A₁ ≈ 1.7 × 10⁻⁴⁴ cm erg⁻² s⁻¹ and A₂ ≈ 1.6 × 10⁻⁵⁶ cm erg⁻² s⁻¹.

By performing an integration over the black hole and the outer region of the accretion disc, the total neutrino annihilation luminosity can be written as

$$\begin{array}{c}\hfill L_{\nu \overline{\nu }}=4\pi {\displaystyle \int \int }{l}_{\nu \overline{\nu }}r\mathrm{d}r\mathrm{d}z.\end{array}$$(19)

This method is physically well founded. However, due to the dependence on the accretion rate, black hole mass and spin parameters, an analytical derivation is not feasible. Consequently, many empirical fitting formulae have been proposed for the neutrino annihilation luminosity. For example, Popham et al. (1999) and Xue et al. (2013) have suggested various parameterised expressions, while some studies have provided corresponding analytical approximations (e.g., Blandford & Znajek 1977). The choice of a specific formula should be based on the range of accretion rates and the research objectives to determine the most suitable fitting function.

The results of Lei et al. (2017) are consistent with the findings of Xue et al. (2013) and Zalamea & Beloborodov (2011) with low and moderate accretion rates. This study adopts the segmented power-law fitting formula proposed by Lei et al. (2017), which better reflects the true trend of annihilation luminosity in the range of Ṁ = 0.01 − 10 M_⊙ s⁻¹, with the black hole mass ranging from 3 M_⊙ to 10 M_⊙. This fitting is particularly suited for simulating the evolution of the central engine energy output during the prompt phase of GRBs.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill L_{\nu \overline{\nu }}& \simeq L_{\nu \overline{\nu },\text{ign}}{[{\left(\frac{\dot{m}}{{\dot{m}}_{\text{ign}}}\right)}^{-{\alpha }_{\nu \overline{\nu }}}+{\left(\frac{\dot{m}}{{\dot{m}}_{\text{ign}}}\right)}^{-{\beta }_{\nu \overline{\nu }}}]}^{-1}\hfill \\ \hfill & \times {[1+{\left(\frac{\dot{m}}{{\dot{m}}_{\text{trap}}}\right)}^{{\beta }_{\nu \overline{\nu }}-{\gamma }_{\nu \overline{\nu }}}]}^{-1},\hfill \end{array}\end{array}$$(20)

where

$$\begin{array}{c}\hfill L_{\nu \overline{\nu },\text{ign}}={10}^{(48.0+0.15a_{\bullet })}{\left(\frac{m_{\bullet }}{3}\right)}^{log\left(\frac{\dot{m}}{{\dot{m}}_{\text{ign}}}\right)-3.3}\text{erg}{\text{s}}^{-1},\end{array}$$(21)

with ${\alpha }_{\nu \overline{\nu }}=4.7$, ${\beta }_{\nu \overline{\nu }}=2.23$, ${\gamma }_{\nu \overline{\nu }}=0.3$, ṁ_ign = 0.07 − 0.063 a_∙, and ${\dot{m}}_{\text{trap}}=6.0-4.0a_{\bullet }^3$.

Similarly, the efficiency of the NDAF and its relation to the luminosity of the observed GRBs can be written as ${\eta }_{\nu \overline{\nu }}{L}_{\nu \overline{\nu }}=f_{\mathrm{b}}{L}_{\gamma ,\text{iso}}$.

It should be noted that the fitting parameters include the luminosity normalisation terms Ė_ν,ign, ${\dot{E}}_{\nu \overline{\nu },\text{ign}}$, the power-law evolution indices α, β, γ, as well as the critical accretion rates ṁ_ign, ṁ_trap. These values are derived from the fits under the conditions of black hole mass M = 3 M_⊙ and spin a_∙ = 0.1 − 0.95.

In the absence of the BZ mechanism, the evolution equation of the Kerr black hole can be derived from the conservation of energy and angular momentum:

$$\begin{array}{c}\hfill \frac{\mathrm{d}{M}_{\bullet }{c}^2}{\mathit{dt}}=\dot{M}{c}^2{E}_{\text{ms}},\end{array}$$(22)

$$\begin{array}{c}\hfill \frac{\mathrm{d}{J}_{\bullet }}{\mathrm{d}t}=\dot{M}{L}_{\text{ms}}.\end{array}$$(23)

Using these physical quantities, the evolution equation of the black hole spin parameter, driven by the neutrino annihilation mechanism, can be derived as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{a}_{\bullet }}{\mathrm{d}t}=\frac{\dot{M}{L}_{\text{ms}}c}{G{M}_{\bullet }^2}-2a_{\bullet }\frac{\dot{M}{E}_{\text{ms}}}{M_{\bullet }}.\end{array}$$(24)

2.3. Comparison of theoretical light curves from two models

To investigate the influence of different central engine mechanisms on the optical behaviour of GRBs during the prompt phase, we adopt the simplified annular accretion disc model proposed by Kumar et al. (2008a); Kumar et al. (2008b), Metzger et al. (2008), Lei et al. (2017), which is regarded as a ring structure concentrated at a representative radius r_d. Under this assumption, the unit mass angular momentum of the disc is calculated using classical orbital mechanics, while the total mass M_d and angular momentum J_d of the disc are equivalently distributed across the ring. The following expression can be derived from the conservation of angular momentum:

$$\begin{array}{c}\hfill j(r_{\mathrm{d}})={\left(G{M}_{\bullet }{r}_{\mathrm{d}}\right)}^{1/2}=\frac{J_{\mathrm{d}}}{M_{\mathrm{d}}}.\end{array}$$(25)

To further obtain the instantaneous accretion rate Ṁ, we relate it to the viscous accretion timescale as

$$\begin{array}{c}\hfill \dot{M}=\frac{M_{\mathrm{d}}}{t_{\text{acc}}}.\end{array}$$(26)

Here, Ṁ denotes the mass accretion rate, and the accretion timescale t_acc is given by $t_{\text{acc}}=\frac{r_{\mathrm{d}}^2}{\nu }\sim \frac{2}{\alpha {\mathrm{\Omega }}_{\mathrm{K}}}$ (Shakura & Sunyaev 1973), where r_d is the radial distance, ν is the kinematic viscosity, and α = 0.1 is the viscosity coefficient. The angular velocity Ω_K, based on the Keplerian velocity under the pseudo-Newtonian potential, is expressed as ${\mathrm{\Omega }}_{\mathrm{K}}={\left(\frac{\mathit{GM}}{R}\right)}^{1/2}\frac{1}{R-R_{\mathrm{g}}}$, where G is the gravitational constant, M is the mass of the central object, R is the radius, and R_g is the Schwarzschild radius.

During the accretion process, the mass and angular momentum of the disc decay and evolve over time, which can be expressed as

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{d}}=-\dot{M},\end{array}$$(27)

$$\begin{array}{c}\hfill {\dot{J}}_{\mathrm{d}}=-L_{\text{ms}}\dot{M}.\end{array}$$(28)

By combining the luminosity functions $L_{\nu \overline{\nu }}(t)$ and L_BZ(t) derived in Sections 2.1 and 2.2, we establish the coupled evolution equations for the black hole and the accretion disc based on the conservation of angular momentum. Solving for the time-dependent evolution of the spin parameter a_∙(t), black hole mass M(t) and accretion rate Ṁ, we can obtain the theoretical evolution of the luminosity light curve driven by the central engine for different initial parameters.

Figures 1 and 2 present the simulated GRB light curves for varying initial parameters, including black hole mass M, spin parameter a_∙, accretion disc mass m_d, accretion rate Ṁ, radiation efficiency ϵ, and beaming factor f_b. The initial conditions are based on different types of compact object mergers or collisions. In merger scenarios, black holes formed from double neutron star mergers typically have masses less than 4 M_⊙, while black holes merging with neutron stars undergo tidal disruption, leading to the formation of an accretion disc, with resulting black hole masses between 4–8 M_⊙ (Shibata et al. 2009). For black holes formed via the collapse of massive stars, the initial mass is typically within the range of 3–10 M_⊙ (Song et al. 2016). In these cases, the initial black hole mass has a minimal effect on GRB production, with the accretion disc mass and accretion process being the main contributors (Lei & Zhang 2011; Qu & Liu 2022). The mass of accretion discs formed in double compact object mergers ranges from 0.1–0.5 M_⊙ (Lee 2001; Fan & Wei 2011; Liu et al. 2012; Lü et al. 2015), while for massive star collapses, it ranges from 1–5 M_⊙ (Zhang et al. 2003; Song et al. 2016). Thus, we set the initial parameter ranges as follows: for long GRBs, m:(3, 10); ṁ:(0.01, 3); m_d:(1, 5); and for short GRBs, m:(2.3, 8); ṁ:(0.01, 3); m_d:(0.1, 0.5). For typical cases, we adopt a radiation efficiency of ϵ = 0.2, a beaming factor of f_b = 0.01, and the black hole spin parameter is taken as a_∙ = (0, 0.5, 0.95). The study of spin equilibrium by Lei et al. (2017) found a_eq ∼ 0.87, while our choice of a_∙ = 0.95 is motivated by a rapidly rotating progenitor, which is necessary for powering a luminous jet at the GRB onset, before the black hole evolves towards equilibrium.

Fig. 1. Simulated long-GRB light curves for two central engine mechanisms: NDAF and BZ. The NDAF model exhibits a steep decay, whereas the BZ model shows a more gradual one. Colours represent variations in initial parameters – (a) black hole spin, (b) accretion rate, (c) accretion disc mass, and (d) black hole mass – illustrated by the colour bars. For each parameter set, the top left panel shows the evolution of the accretion rate over time; the top right panel shows the evolution of the spin parameter; the bottom left panel shows the luminosity evolution, where the NDAF model is plotted as red dashed lines (typical slope −4.5) and the BZ model as blue dashed lines (typical slope −1.5); and the bottom right panel compares the two mechanisms after alignment.

Fig. 2. Same as Fig. 1 but for simulated short-GRB light curves.

By fitting the decay slope of the theoretical light curve, we find that on the logarithmic time scale, the light curves dominated by the two mechanisms exhibit distinct differences in the slope of the decay phase for both long and short GRBs. The details are as follows:

1. When the luminosity varies by several orders of magnitude, the slope of the decay phase effectively differentiates between the two mechanisms.

2. A steeper luminosity decay in the NDAF mechanism, with slopes from −4.05 ± 0.08 to −7.01 ± 0.40, where the typical luminosity decays as L ∝ t^−4.5. Theoretically, the expected evolution of the NDAF luminosity depends on the accretion regime. For high accretion rates (ṁ ≫ ṁ_ign), the luminosity scales as $L_{\nu \overline{\nu }}\propto {\dot{M}}^{2.23}\propto t^{-3.72}$. For low accretion rates (ṁ ≪ ṁ_ign), it scales as $L_{\nu \overline{\nu }}\propto {\dot{M}}^{4.7}\propto t^{-7.83}$. The decay slopes range approximately from t^−3.7 to t^−7.8, consistent with the fitted slopes.

3. The luminosity decay for the BZ mechanism is relatively shallow, with slopes typically ranging from −1.08 ± 0.04 to −1.49 ± 0.05, consistent with L ∝ t^−1.5. Theoretically, the BZ luminosity is linearly proportional to the accretion rate, thus predicting a decay slope of L_BZ ∝ Ṁ ∝ t^−5/3 ≈ t^−1.67, closely matching the fitted results.

This phenomenon suggests that the slope of the light curve during the decay phase serves as a critical parameter for distinguishing between the different mechanisms of the central engine of GRBs. We next integrate actual data from single-pulse GRBs exhibiting FRED structures with the theoretical curves outlined above to perform a fitting analysis, with the aim of evaluating the model’s potential for differentiating these mechanisms.

3. Samples and methods

To ensure clear and well-defined light curves with reliable fitting, this study focused on single-peaked GRB events that exhibit a typical fast rise exponential decay (FRED) profile. Here, a “single-peaked GRB” refers to the emergent overall profile that results from the superposition of multiple local dissipation episodes throughout the GRB duration, with each episode contributing to the overall pulse. We assume that the peak flux of this variation is proportional to the instantaneous energy release of the central engine. In this case, the overall pulse profile primarily traces the temporal evolution of the central engine. This provides a physical basis for our statistical study of the GRB central engine based on the pulse profile. This simple structure improves fitting accuracy and facilitates mechanism identification, while also aiding the identification and interpretation of the underlying physical mechanisms.

To connect jet power with observed pulse shapes, we derive the observed variability timescale t_obs based on radius R, Lorentz factor Γ, and redshift z of the emitting region. The relationship is given by $t_{\text{obs}}=\frac{R(1+z)}{2c{\mathrm{\Gamma }}^2}$. For typical GRB parameters, such as Γ ∼ 100−300 and R ∼ 10¹³ cm, the timescale is approximately 0.01 s. If the radiation region extends to R ∼ 10¹⁴ cm, t_obs would be about 0.17 s. This suggests that even a single-peaked pulse contains multiple dissipation events, with the overall decay trend linked to energy injection from the central engine. Thus, the decay provides a reliable indicator of the central engine’s temporal power evolution.

The light curve data used in this study are from the Burst Alert Telescope (BAT) on the Swift satellite from December 2004 to May 2025 (Lien et al. 2016). This instrument provides multi-channel count rate data that have been processed with mask-weighting, offering high time resolution. The four energy bands of BAT are as follows: 15–25 keV, 25–50 keV, 50–100 keV, and 100–350 keV. To reflect the total luminosity variation of the GRB, we combined the data from all four energy bands and analysed the broadband light curves. Additionally, Swift offers light curve data with varying time resolutions. By comparing these light curves, we found that the 2 ms and 8 ms resolution light curves are suitable for short GRBs. However, the 1 s resolution light curves, with fewer data points in the pulse regions, result in rougher curves and may exhibit complex pulse overlaps. Therefore, we selected 64 ms resolution data to identify and analyse the FRED single-pulse structure.

We selected GRB events with single-pulse structures and 64 ms resolution from a specific database¹. A total of 85 single-pulse samples were collected, including 39 with redshift. For GRBs lacking redshift, we estimate the pseudo-redshift using the standard ΛCDM cosmological parameters (H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.27, Ω_Λ = 0.73 (Spergel et al. 2003)) in combination with the Yonetoku relation. This approach ensures that the calculated luminosity evolution of the light curves can be compared with the theoretical light curves.

To analyse the structural characteristics of GRB light curves in the prompt phase, we used the MCMC method with the empirical Kocevski–Ryde–Liang (KRL) function to fit the pulses and evaluate the goodness of fit for the FRED profile using the minimum chi-square method (Kocevski et al. 2003; Liang et al. 2010; Foreman-Mackey et al. 2013), which is defined as

$$\begin{array}{c}\hfill F(t)=F_{\mathrm{p}}{\left(\frac{t+t_0}{t_{\mathrm{p}}+t_0}\right)}^r{[\frac{d}{r+d}+\frac{r}{r+d}{\left(\frac{t+t_0}{t_{\mathrm{p}}+t_0}\right)}^{r+1}]}^{-\frac{r+d}{r+1}}.\end{array}$$(29)

The KRL function consists of five adjustable parameters that allow for flexible fitting of light curves with various shapes. Specifically, F_p represents the peak luminosity, t_p is the peak time of the pulse, and r (d) are the power-law indices for the rise (decay) phase, respectively. This equation is valid for t > −t₀, where t₀ represents the offset of the pulse start time relative to the trigger time. The temporal boundaries of the pulse (t_start and t_end) are determined by the signal-to-noise ratio threshold (S/N = 1), corresponding to the initial and final epochs when the pulse profile intersects the detection threshold.

To determine t₀, we adopted an empirical method similar to that of Kocevski et al. (2003), in which 10% of the standard deviation of the background noise is added to its mean, i.e., t₀ = μ+ 0.1 × σ, where μ is the mean and σ is the standard deviation of the background noise. The signal onset is defined as the point where the signal first exceeds this threshold. This method simplifies the KRL function to four adjustable parameters. By applying the t₀-corrected KRL function, we can directly compare its parameters with those of theoretical models, such as the BZ mechanism or the NDAF models, which assume that the light curve evolves from t = 0 without accounting for the trigger delay.

The data processing procedure was as follows:

1. The count rates (counts s⁻¹ det⁻¹)² for each energy channel are converted to energy flux (erg cm⁻² s⁻¹) using the best-fitting model (power law or broken power law).

2. For samples lacking redshift, we estimate the pseudo-redshift using the Yonetoku relation. The flux values from the four energy bands are then combined, and K-correction is applied to obtain the light curve showing the luminosity evolution over time.

3. We determine the pulse start time for each sample and apply MCMC with the KRL model (with t₀ fixed) to fit the pulses. The observed time T_obs is converted to rest-frame time $T_{\mathrm{rest}}=\frac{T_{\mathrm{obs}}}{1+z}$ (with z being the source redshift), and the decay slope parameter d is derived.

The spectral fitting in this study used both the power law (PL) and broken power law (CPL) models. The specific forms of these models are as follows:

$$\begin{array}{c}\hfill N_{E,\mathrm{CPL}}(E)=N_{0,\mathrm{CPL}}{\left(\frac{E}{100\mathrm{keV}}\right)}^{{\alpha }_{\mathrm{CPL}}}exp(-\frac{E}{E_{\mathrm{p}}}),\end{array}$$(30)

$$\begin{array}{c}\hfill N_{E,\mathrm{PL}}(E)=N_{0,\mathrm{PL}}{\left(\frac{E}{100\mathrm{keV}}\right)}^{{\alpha }_{\mathrm{PL}}},\end{array}$$(31)

where N_ECPL(E) and N_E,PL(E) represent the differential photon spectra for the CPL and PL models, respectively (Schaefer 2007). Here, α_CPL and α_PL are the spectral indices, and E_p is the break energy for the CPL model.

The count rate is converted into flux using the following equation:

$$\begin{array}{c}\hfill F_{\mathrm{p}}=P_{\mathrm{p}}\frac{{\int }_{E_{\text{min}}}^{E_{\text{max}}}EN(E)\mathrm{d}E}{{\int }_{E_{\text{min}}}^{E_{\text{max}}}N(E)\mathrm{d}E},\end{array}$$(32)

where F_p is the flux (erg cm⁻² s⁻¹), P_p is the photon peak flux (photons cm⁻² s⁻¹), and the integrals are performed over the energy range [E_min, E_max]. These integrals are evaluated over the energy intervals 15–25 keV, 25–50 keV, 50–100 keV, and 100–350 keV. The total flux is then computed by summing the contributions from each range.

The isotropic luminosity L_iso is then calculated as

$$\begin{array}{c}\hfill L_{\text{iso}}=4\pi D_{\mathrm{L}}^2{F}_{\mathrm{p}}k,\end{array}$$(33)

where D_L is the luminosity distance and k is the K-correction factor, which accounts for the cosmological redshift and is given by

$$\begin{array}{c}\hfill k=\frac{{\int }_{1/(1+z)}^{{10}^4/(1+z)}EN(E)\mathrm{d}E}{{\int }_{E_{\text{min}}}^{E_{\text{max}}}EN(E)\mathrm{d}E},\end{array}$$(34)

where z is the redshift, [E_min, E_max] = 15–350 keV, and both integrals are evaluated in the source frame.

We estimate the spectral peak energy E_p using the empirical E_p − L correlation (Yonetoku et al. 2004; Nava et al. 2012)

$$\begin{array}{c}\hfill log\left[E_{\mathrm{p}}(1+z)\right]=-25.33+0.53logL,\end{array}$$(35)

where L is the isotropic luminosity and z is the redshift.

4. Results

To investigate the differences in light-curve morphology between NDAF and BZ mechanisms, we performed spectral fitting using the KRL function and analysed 85 single-pulse GRBs. The results of the light-curve fits, along with the isotropic gamma-ray energy (E_γ,iso) and duration (T₉₀), are listed in Table 1 (available in machine-readable form).

Figure 3 compares the observed data (black) with the fitted KRL models (red lines), highlighting the decay behaviours characteristic of each mechanism. The dotted lines indicate the pulse start and end times, and the data have been corrected for redshift. The energy range for the fitting was initially 15–350 keV and was then extended to the full energy range of 1–10 000 keV through K-correction. Each panel represents a unique GRB, illustrating the applicability of the KRL model across varying attenuation slopes and mechanisms, and providing insights into the physical processes governing the light-curve evolution.

Fig. 3. Examples of fitting light curves (1–10 000 keV) to the KRL model (red lines) for GRB 081222, GRB 161218A, GRB 170101A, GRB 180728A, GRB 200922A, GRB 201105A, GRB 210410A, GRB 211225B, and GRB 250605A. The dashed vertical lines indicate the pulse start and end times. The observed data (shown in grey) have been corrected for redshift, and the corresponding fitting parameters are listed in Table 1.

To further assess the accuracy of the model, Figure 4 compares the theoretical light curves of typical samples with the observed data after aligning their peak values. The x-axis represents the observation time corrected for redshift, i.e., T/(1+z). In the upper panel of Figure 4, the red curve represents GRB 120326A, which exhibits an observed decay slope of d ≈ 4.49, matching the steep luminosity decline predicted for the NDAF mechanism (d ≈ 3.7–7.8). Conversely, the lower panel shows GRB 160131A, with a shallower decay slope of d ≈ 1.59, in close agreement with the BZ mechanism prediction (d ≈ 1.67). Theoretical and observational results are in strong agreement, providing robust evidence for these mechanisms in explaining the observed GRB light curves.

Fig. 4. Comparison between theoretical models and observational data for (a) GRB 120326A and (b) GRB 160131A. Theoretical light curves are shown for the Blandford–Znajek (BZ) mechanism (solid lines) and the neutrino-dominated accretion flow (NDAF) mechanism (dashed lines), following the same convention as in Figures 1 and 2. Observational data are indicated by red circles, with lighter shading showing the error bars. The upper panel exhibits a steep luminosity decay ($L_{\nu \overline{\nu }}\propto t^{-3.72}$–t−7.83), consistent with NDAF predictions, whereas the lower panel shows a shallower decay (L ∝ t−1.67), matching BZ model expectations.

In a sample of 85 GRBs, we find that the decay slopes of 15 bursts (e.g., GRB 050717, GRB 070318, GRB 070808, GRB 080805, and GRB 090530) are consistent with the BZ mechanism, while 22 bursts (e.g., GRB 081222, GRB 090129, GRB 091018, GRB 110318A, GRB 120213A, and GRB 120326A) are more in line with the NDAF mechanism (see Table 1). However, based on the distribution of decay slopes across all samples (Figure 5), most bursts cluster in the range 2 < d < 4. This range lies between the typical values predicted by the NDAF (d ≈ 3.7–7.8) and BZ (d ≈ 1.67) mechanisms shown in Figures 1 and 2, indicating that many events cannot be classified into a single mechanism solely from their decay slopes. For the statistical analysis in Figure 5, we excluded three extreme cases (GRB 070306, GRB 140209A, and GRB 161004B) with d > 7.8, as such steep slopes are likely influenced by observational uncertainties or atypical physical conditions, and would disproportionately skew the high-d tail of the distribution. A Gaussian mixture model (GMM) fit to the distribution yields three peaks at d ≈ 2.04 (weight 0.50), d ≈ 3.56 (weight 0.32) and d ≈ 6.04 (weight 0.18). The first two peaks fall within 2 < d < 4, while the third represents the high-d tail of the distribution.

Fig. 5. Distribution of decay slopes d for 85 GRBs, fitted with a three-component Gaussian mixture model (GMM). Three extreme cases (GRB 070306, GRB 140209A and GRB 161004B) with d > 7.8 were excluded to avoid skewing the high-d tail due to potential observational or physical anomalies, and their removal does not affect the main results. The best-fit components peak at d ≈ 2.04, d ≈ 3.56 and d ≈ 6.04 with relative weights of 0.50, 0.32 and 0.18, respectively. Most events fall in the range 2 < d < 4, suggesting hybrid central engines, differing from idealised BZ (d ≈ 1.67) or NDAF (d ≈ 3.7–7.8) predictions.

5. Discussion

The interpretation of GRB pulse profiles should account for both local and global physical mechanisms. Models involving high-latitude emission from a spherical shell or localized magnetic reconnection events (Genet & Granot 2009; Uhm & Zhang 2015, 2016; Uhm et al. 2024) can produce a decaying trend without invoking a decaying jet power, and these processes likely govern the fine-scale temporal structure (e.g., micro-spikes) observed within pulses. However, in this work, we focus on the macroscopic envelope of the pulse, which is primarily influenced by the global energy injection history of the central engine. We suggest that the global decay trend, particularly over timescales longer than those governed by local dissipation processes, may trace the temporal evolution of the jet power. While local processes, such as high-latitude emission and magnetic reconnection, govern the rapid fluctuations within the pulse, the overall decay remains coupled to the global energy injection history of the central engine. We argue that the overall decay trend remains statistically indicative of the temporal evolution of the central engine.

This framework motivates the analysis of the temporal decay slope as a proxy for the engine’s behavior. The observed concentration of intermediate slopes (2 < d < 4) may reflect hybrid central engine mechanisms or transitional accretion regimes operating in a significant fraction of bursts. Assuming that the evolution of the light curve is primarily driven by accretion processes, this slope can be used to differentiate between the two mechanisms (see Figures 1 and 2). When the photometric attenuation does not reach a full order of magnitude, the measured slope is often intermediate between the two canonical values. In such cases, the trend in our sample leans toward NDAF-like behaviour. Instrumental factors such as orientation effects or calibration uncertainties can also introduce data gaps that obscure the intrinsic decay pattern. For bright, long-duration events with high-quality light curves, however, the slope remains a useful diagnostic for distinguishing between the two mechanisms.

A Gaussian mixture model (GMM) fit to the slope distribution supports this interpretation: while a two-component model is favoured by the Bayesian information criterion (BIC), a three-component fit reveals peaks near d ≈ 2.04, 3.56 and 6.40, with the high-d tail likely tracing rare NDAF-dominated bursts or specific geometric and radiative configurations (Figure 5). The concentration of slopes in the 2 < d < 4 range persists regardless of whether two- or three-component fits are adopted, reinforcing the prevalence of intermediate decay behaviours. Physical processes beyond the idealised models may further shape the observed slopes. Strong disc outflows, for instance, can remove a substantial fraction of the accreting mass before it reaches the black hole (Liu et al. 2018), lowering the late-time accretion rate and steepening the decay.

The light-curve slope may be influenced by magnetic field variations during the BZ process, such as attenuation or sudden disappearance of the magnetic field, which could increase the slope steepness. While magnetic flux is assumed constant for simplicity, variations may affect the slope. We identify GRBs such as GRB 170317A, GRB 141004A, GRB 120521C, GRB 120326A, GRB 100814A, GRB 100704A, GRB 080430, GRB 071003, GRB 070306, GRB 061202, and GRB 060904B, which exhibit plateaus and a steep decay for afterglows, indicative of the magnetar model (Dai 2004; Zhang et al. 2006). Magnetic field evolution could alter the prompt light-curve structure, influencing intensity and spectral features. The correlation between the prompt and afterglow phases requires further investigation through multi-band data (optical, X-ray, radio) (Yi et al. 2016, 2021, 2022) and numerical simulations. Our analysis focuses on the prompt phase; if the afterglow phase involves continued accretion or fallback, the actual accreted mass could be higher, consistent with conclusions drawn from X-ray afterglow analyses (e.g., Wu et al. 2013; Hou et al. 2014; Gao et al. 2016). Additionally, the complexity of GRB pulse structures should be considered. Some GRB pulses may result from the superposition of multiple components, complicating light-curve analysis and making mechanism differentiation more difficult. Small peak structures, typically generated by internal shock wave collisions, can subtly influence the overall light-curve behaviour, despite being unrelated to the primary interaction. The choice of binning, fitting interval range and initial time point t₀ may also influence the slope calculation, though their effects are generally minimal.

In summary, our study provides new insights into GRB central engine mechanisms through light-curve slope analysis. Under low photometric attenuation, slopes in the transition zone between the BZ and NDAF mechanisms are harder to distinguish, yet well-sampled, bright events retain clear diagnostic value. Targeted multi-wavelength campaigns, combined with prompt–afterglow correlations, will be essential to connect slope evolution with magnetic field dynamics and the underlying engine activity.

6. Conclusions

We propose a method to distinguish between the BZ and NDAF mechanisms through analysis of GRB light-curve decay slopes. Our results show that the BZ mechanism typically follows L_BZ ∝ t^−1.67, while the NDAF mechanism exhibits slopes ranging from $L_{\nu \overline{\nu }}\propto t^{-3.7}$ to t^−7.8, highlighting the decay slope as a key diagnostic. From 85 Swift/FRED single-pulse GRB samples at 64 ms resolution, 15 events are consistent with the BZ mechanism and 22 with the NDAF mechanism. These results suggest that the NDAF mechanism explains the majority of GRBs in our sample, complementing previous studies that indicate high-luminosity GRBs may require the BZ mechanism (Liu et al. 2015; Yi et al. 2017). The contrasting decay behaviours align with differences in accretion processes: the BZ mechanism produces smoother decays via gradual energy extraction, whereas the NDAF mechanism yields steeper decays due to rapid accretion rate decline in the present model framework.

While these findings establish a potential link between the decay slope and central engine mechanism, they also motivate further investigation. Several limitations remain, including the thin-disc approximation and the assumption of purely single-mechanism models. Future studies should relax these assumptions and incorporate multi-band observations, together with correlations between X-ray afterglows, gravitational wave signals, and central engine activity. Model refinements that account for hybrid mechanisms, disc precession and magnetic reconnection, supported by new telescopes, improved calibration and multi-band data, combined with numerical simulations, will be essential for a more complete understanding of GRB light-curve formation, enhancing the accuracy of mechanism identification and deepening our understanding of the underlying GRB physics.

Data availability

Table 1 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/707/A392.

Acknowledgments

We thank Xiu-Juan Li, Wen-Long Zhang and Cui-Ying Song for the helpful discussions. This work is supported by Shandong Provincial Natural Science Foundation (ZR2025MS47, ZR2021MA021), the National Natural Science Foundation of China (Grant Nos. 12494575, 12473012 and 12494572), the Natural Science Foundation of Jiangxi Province of China (grant No. 20242BAB26012), the National Natural Science Foundation of China under grants 12473012 and 12533005, and the National Key R&D Program of China (No. 2023YFC2205901) and Manned Spaced Project (CMS-CSST-2021-A12).

References

Appendix A: Figures

Fig. A.1. Light-curve fitting of 85 GRBs using the KRL function (red lines). The dotted lines indicate the pulse start and end times. The data have been corrected for redshift, and the observed data are shown in grey.

All Figures

Fig. 1. Simulated long-GRB light curves for two central engine mechanisms: NDAF and BZ. The NDAF model exhibits a steep decay, whereas the BZ model shows a more gradual one. Colours represent variations in initial parameters – (a) black hole spin, (b) accretion rate, (c) accretion disc mass, and (d) black hole mass – illustrated by the colour bars. For each parameter set, the top left panel shows the evolution of the accretion rate over time; the top right panel shows the evolution of the spin parameter; the bottom left panel shows the luminosity evolution, where the NDAF model is plotted as red dashed lines (typical slope −4.5) and the BZ model as blue dashed lines (typical slope −1.5); and the bottom right panel compares the two mechanisms after alignment.

In the text

	Fig. 2. Same as Fig. 1 but for simulated short-GRB light curves.
In the text

Fig. 3. Examples of fitting light curves (1–10 000 keV) to the KRL model (red lines) for GRB 081222, GRB 161218A, GRB 170101A, GRB 180728A, GRB 200922A, GRB 201105A, GRB 210410A, GRB 211225B, and GRB 250605A. The dashed vertical lines indicate the pulse start and end times. The observed data (shown in grey) have been corrected for redshift, and the corresponding fitting parameters are listed in Table 1.

In the text

Fig. 4. Comparison between theoretical models and observational data for (a) GRB 120326A and (b) GRB 160131A. Theoretical light curves are shown for the Blandford–Znajek (BZ) mechanism (solid lines) and the neutrino-dominated accretion flow (NDAF) mechanism (dashed lines), following the same convention as in Figures 1 and 2. Observational data are indicated by red circles, with lighter shading showing the error bars. The upper panel exhibits a steep luminosity decay ($L_{\nu \overline{\nu }}\propto t^{-3.72}$–t−7.83), consistent with NDAF predictions, whereas the lower panel shows a shallower decay (L ∝ t−1.67), matching BZ model expectations.

In the text

Fig. 5. Distribution of decay slopes d for 85 GRBs, fitted with a three-component Gaussian mixture model (GMM). Three extreme cases (GRB 070306, GRB 140209A and GRB 161004B) with d > 7.8 were excluded to avoid skewing the high-d tail due to potential observational or physical anomalies, and their removal does not affect the main results. The best-fit components peak at d ≈ 2.04, d ≈ 3.56 and d ≈ 6.04 with relative weights of 0.50, 0.32 and 0.18, respectively. Most events fall in the range 2 < d < 4, suggesting hybrid central engines, differing from idealised BZ (d ≈ 1.67) or NDAF (d ≈ 3.7–7.8) predictions.

In the text

In the text

	Fig. A.1. Light-curve fitting of 85 GRBs using the KRL function (red lines). The dotted lines indicate the pulse start and end times. The data have been corrected for redshift, and the observed data are shown in grey.
In the text