© The Authors 2023

1. Introduction

Gamma-ray bursts (GRBs) are among the most fascinating transients in astrophysics, and the recent association of the gravitational wave source GW170817 with the short GRB 170817A (Abbott et al. 2017a,b,c; Goldstein et al. 2017; Savchenko et al. 2017) has further increased the interest in these objects. GRBs are hence not only a unique laboratory for high-energy astrophysical phenomena, but also an ideal multi-messenger factory. Although enormous advances have been made over the past two decades in understanding their emission mechanisms, many questions still need to be answered.

The fast re-pointing capabilities (on time scales of minutes) of the X-Ray Telescope (XRT) on board the Neil Gehrels Swift Observatory (Swift hereafter, Gehrels et al. 2004) has enabled the discovery of several features not predicted by the standard afterglow model (e.g. Meszaros & Rees 1993b,a; Sari et al. 1998) in the early GRB afterglow light curve stages (< 0.5 day). Before Swift, GRB afterglows were known as sources of non-thermal radiation quickly fading with time as power laws (F ∝ t^−a with a > 0.8 − 1). Today, it is known that soon after the end of the prompt emission, a large fraction of GRB X-ray afterglow light curves show a peculiar shallow flux decay phase (‘plateau’; e.g. Nousek et al. 2006; Zhang et al. 2006; O’Brien et al. 2006) during which the temporal decay index reaches values in the range a ∈ [0, ∼0.7] for a temporal interval that typically lasts for a few hours. Post-plateau fluxes decay according to standard afterglow predictions, and the end of the plateau phase is typically not accompanied by any spectral evolution. These properties are not compatible with the standard afterglow model based on synchrotron emission from a shock-accelerated electron population (e.g. Sari et al. 1998; Granot & Sari 2002). So far, no firm conclusion has been reached on the origin of plateaus, but it is commonly believed that they encode crucial information on the GRB central engine.

The leading interpretation proposed for the origin of plateaus invokes the presence of a central engine that can provide continuous energy injection, for example, fallback onto a newly formed black hole or dipole radiation from magnetar spin down. Both of these scenarios can fairly well describe the observed light curve morphology (e.g. Rowlinson et al. 2013; Li et al. 2018). An important step forward was recently made by Stratta et al. (2018). They successfully compared the magnetar model with a data set of 51 GRBs, finding evidence of consistency of the GRB distribution in the diagram of magnetic field versus spin period (B–P) with the well known spin-up line for accreting Galactic X-ray pulsars (e.g. Bhattacharya & van den Heuvel 1991; Pan et al. 2013). In particular, they report that the normalisation of the observed relation perfectly matches the spin-up line predictions for typical neutron star masses (∼1 M_⊙) and radii (∼10 km) and for mass accretion rates expected in GRBs (10⁻⁴ M_⊙ s⁻¹ < Ṁ < 0.1 M_⊙ s⁻¹). This result was independently confirmed two years later (Lin et al. 2020), supporting the presence of a magnetar as the central engine.

Another interpretation invokes a completely different scenario where the plateau is created by the high latitude emission (HLE; Fenimore et al. 1996) from a structured jet (SJ) (Dai & Gou 2001; Lipunov et al. 2001; Rossi et al. 2002) associated with the prompt phase of the GRB (Oganesyan et al. 2020). The HLE consists in the radiation observed from larger angles relative to the line of sight after the prompt emission from a curved surface is switched off (Kumar & Panaitescu 2000). Before GW170817/GRB 170817A, the GRB jet structure was typically approximated with a uniform energy and velocity distribution within the cone aperture. The off-axis view of GRB 170817A has enabled us for the first time to observe evidence of the structured nature of the ultra-relativistic jet, where energy and velocity distribution decrease gradually toward the jet edges. Independent of the jet structure, the HLE is typically invoked to explain the initial steep decay phase of X-ray light curves (e.g. Liang et al. 2006; Tagliaferri et al. 2005; Barthelmy et al. 2005) or X-ray flares (e.g. Rossi et al. 2011). In the case of an SJ, the expected HLE lasts longer. Indeed, compared to a top-hat jet, photons departing from a given angular distance from the jet axis (θ) are less Doppler boosted, and therefore the radiation is less beamed. This effect gives an extra contribution to the observed flux at late times, explaining the flattening and the longer duration of the HLE. This results in a shallow segment of the light curve for a portion of time consistent with the duration of the observed plateau. The HLE from an SJ has the ability to reproduce the temporal evolution of the plateau as well as its spectral properties.

Beniamini et al. (2020a) recently proposed a different interpretation where X-ray plateaus can form during the afterglow phase if the observer is slightly off axis from the jet, though no direct comparison with observed light curves are performed. They show how this interpretation can explain most of the main plateau properties, including a correlation between the prompt energy, the plateau luminosity, and duration, which is specific to their model and not expected in energy injection scenarios. Since all these interpretations can fairly well reproduce the plateau features observed in X-rays, it is challenging to identify the actual mechanism.

This work is the continuation of a previous work (Stratta et al. 2022) where we explored the behaviour of the plateaus over a larger energy band, in particular adding optical observations, to investigate if a specific scenario is prefererred. A considerable fraction of events also show a plateau feature at optical wavelengths (e.g. Si et al. 2018; Dainotti et al. 2020), but a comprehensive consistency check of X-ray and optical light curves within the afterglow model is still confined to a few cases (e.g. Gompertz et al. 2015; Zhang et al. 2018). In this work, we use a sample of 30 GRBs with evidence of an X-ray plateau and simultaneous optical data during and after the plateau phase in order to test different scenarios, including the energy injection from a magnetar, the HLE from an SJ and modifications to the standard forward shock model. In Sect. 2, we present the sample and in Sect. 3, we present our analysis of the X-rays and optical plateaus. We discuss our results in Sect. 4 and draw our conclusions in Sect. 5.

2. GRB sample

2.1. X-ray plateau sample selection

We considered all the GRBs with X-ray afterglow detected with Swift from 2005 up to mid 2018. In order to identify GRBs with evidence of a plateau feature in their afterglow, we used the publicly available Swift XRT Repository¹ (Evans et al. 2007, 2009) and the provided data analysis tools to fit multiple power laws along the light curve of each GRB. We selected all the GRBs that present at least one segment with a temporal slope (−0.8 ≤ α ≤ 0.8 within errors) in the X-ray light curve. This phenomenology is indeed challenging for the standard afterglow model that assumes a constant density environment, which is typically observed in GRBs (e.g. Sari et al. 1998). A wind environment can give rise to a flat temporal evolution, but it requires one of the following conditions: (1) If the jet is already in the deceleration phase, there must be a rising spectral slope in the F_ν − ν plane, which is a condition never observed in X-ray afterglows. (2) The spectral slope can be negative, but the jet has to be in the pre-deceleration phase. The second case would require rather low values of the bulk Lorentz factor (Γ ∼ 5 − 10) since the X-ray plateau can last up to 10⁴ − 10⁵ s. Defining t_i and t_f as the initial and final time of the plateau, respectively, we also discarded all GRBs where t_f/t_i is less than two. This requirement was necessary in order to select a long enough plateau that would have a sufficient number of photons to perform a time-resolved spectral analysis, as explained in Sect. 3.1.

2.2. X-ray plateau optical counterpart sample

In order to investigate on the broad-band properties of the plateau, in this work we focus on a sample subset with multi-band optical-NIR follow-up during the X-ray plateau phase. Specifically, we selected GRBs for which multi-band data during the X-ray plateau and post-plateau phase are enough to perform a time-resolved spectral analysis (see next section).

The optical sample is based on criteria and methods originally presented in Kann et al. (2006). The sample includes GRBs with afterglow, redshift, and good coverage in the UV-optical-NIR domain and is based on the updated sample and analysis of works in preparation (Kann et al., in prep.). For this work in particular, we re-analysed GRB 110213A, while the analysis of GRB 180728A is taken from Rossi et al. (in prep.). Hereby, we assumed achromaticity and that all bands follow the same temporal evolution. Any contribution to the afterglow flux from the presence of the GRB host galaxy and/or an associated supernova was subtracted. The spectral energy distributions (SEDs) were fitted with local dust extinction laws (Pei 1992), determining the intrinsic spectral slope β₀ and the extinction in the rest frame V band A_V (in mag). Using the method of Kann et al. (2006), the light curves were then corrected for line-of-sight extinction. The achromatic nature of the light curves allowed us to shift data from other bands to the common R_C band, creating a maximally dense light curve for each analysed GRB afterglow. We emphasise that we used results based only on the optical-NIR SED. We report the A_V value of each GRB of our sample in Table A.2. Most of these measurements are found in the work presented in Kann et al. (2010, 2011) and Kann et al. (in prep). There are a few cases (indicated in the table) where we referred to works dedicated to single events. Our final sample of GRBs with X-ray plateau and optical counterparts consists of 30 events. General information of the full GRB sample, including burst duration, redshift, and energetics, is presented in Table A.1.

Throughout this work, the flux density of the afterglow is described as F_ν(t)∝t^−αν^−β. We assumed a ΛCDM cosmological model for calculations. All errors are at 1σ unless otherwise specified.

3. Data analysis

Our goal is to verify that both the optical and X-ray emission during the X-ray plateau are consistent with a single synchrotron component in each temporal bin. To achieve this goal, we performed a time-resolved spectral analysis in the XRT band, fitting the X-ray data with a single absorbed power law. From this analysis, we derived the X-ray flux and the photon index. We then extrapolated the expected flux in the optical regime at the X-ray binning times (see Sect. 3.2 and Fig. 1), and we checked whether the extrapolation is compatible with the optical observations shifted to the R_C band (see Sect. 2.2).

Fig. 1. Schematic representation of the method used to compare X-ray and optical data for each temporal bin. The minimum and maximum optical fluxes allowed by the standard afterglow model assuming that there is a single break νopt < νb < νX are indicated by Fm,opt and FM,opt, respectively. We derived FX, Fopt, and the X-ray spectral slope from spectral analysis. Therefore, if Fm,opt < Fopt < FM,opt, we then considered that the optical counterpart is consistent with a single synchrotron spectrum, and we could derive νb. The blue line has no break, the orange line has νb = νX, and the green line has a break between the optical and X-ray bands. We note that our assumption here is that νm < νopt, and it is supported by the observation of βopt > 0 for all the analysed events (see Sect. 3.2) that exclude the case νb = νm.

3.1. X-ray time-resolved spectral analysis

The X-ray light curves were re-binned imposing that the number of counts per bin in the band (0.5–10) keV were to be above a certain threshold N₀, and we chose N₀ > 500 counts. For bins with N₀ < 500, the spectral analysis either gave errors that were too large on the parameters or the fit did not converge at all, due to the noisiness of the spectrum. Since the X-ray light curves contain some observational gaps, the temporal length was not the same for all the bins. Moreover, given the criteria described above, the bin lengths tended to be larger in the post-plateau phase, where the count rate decreases.

We performed a spectral analysis for each bin using XSPEC, version 12.10.1, and PyXspec. We considered only the photons in the band E = 0.5 − 10 keV. Each spectrum was modelled by adopting an absorbed power law, and for the absorption, we used the Tübingen-Boulder model (Wilms et al. 2000). Since the redshift is known for all the GRBs of our sample, we were able to distinguish a Galactic absorber and the host galaxy absorber. The Galactic absorption is taken from Kalberla et al. (2005). The specific syntax in XSPEC was ‘tbabs*ztbabs*po’. The estimation of the host equivalent hydrogen column density N_H was performed on the time-integrated spectrum of the post-plateau phase, where we did not expect strong spectral evolution, as verified by Butler & Kocevski (2007) and Mu et al. (2016). Once the host N_H was determined, it was fixed to be the same for all the bins of the light curve. The only free parameters were the photon index and normalisation. Such a procedure is preferred over the approach of leaving the host N_H as a free parameter because of the degeneracy between the photon index and the column density.

3.2. Comparison of optical and X-ray data

In this section, we show a method to derive the overall spectral properties of the plateau and post-plateau phases, exploiting the knowledge of the X-ray flux, the optical flux, and the X-ray photon index. In the synchrotron scenario, two characteristic frequencies are defined: (1) the cooling frequency ν_c, which is associated with the cooling Lorentz factor ${\gamma }_c=\frac{6\pi m_{\mathrm{e}}c}{{\sigma }_{\mathrm{T}}t{B}^2}$, where the time t and the magnetic field B are defined in the comoving frame, and (2) the frequency ν_m, which is associated with the minimum Lorentz factor γ_m of the electron energy distribution.

In the following, we work in the assumption of a slow cooling regime, namely, ν_m < ν_c, as it is commonly found in GRB afterglows. Indeed, when computing the ratio between the cooling frequency ν_c and ν_m (Sari et al. 1998), we obtain:

$$\begin{array}{c}\hfill \frac{{\nu }_c}{{\nu }_m}\sim {10}^5{\epsilon }_{\mathrm{B},-3}^{-2}{E}_{52}^{-1}{n}^{-1}{\epsilon }_{\mathrm{e},-1}^{-2}{t}_d{(1+Y)}^{-2},\end{array}$$(1)

where ϵ_B and ϵ_e are the fraction of energy that go to the magnetic field and electrons, respectively; E is the isotropic kinetic energy; n is the circumburst density; t_d is the time measured in days; and the Y parameter indicates the relevance of the Synchrotron Self Compton (SSC) component (we adopted the notation Q_X = Q/10^X). In Appendix A, we show how the ratio ν_c/ν_m depends on the Y parameter. We found that for typical micro-physical values, it is more plausible to have ν_c > ν_m, even if the opposite regime cannot be rejected a priori (discussed later in this section). In the slow cooling regime, we expected the flux density to be:

$$\begin{array}{c}\hfill F_{\nu }\propto {\nu }^{-\beta }\end{array}$$

for ν > ν_c and

$$\begin{array}{c}\hfill F_{\nu }\propto {\nu }^{-(\beta -1/2)}\end{array}$$

for ν < ν_c, where β = p/2 and p is the power law index of the particle distribution. If ν_c < ν_opt, then the optical and X-ray fluxes can be connected with a single power law with spectral index β. In this case we define the maximum allowed flux²:

$$\begin{array}{c}\hfill F_{M,\mathrm{opt}}=F_X{({\nu }_{\mathrm{opt}}/{\nu }_X)}^{-\beta }.\end{array}$$

If instead ν_opt < ν_c < ν_X, the optical flux is below the spectral extrapolation from X-rays assuming a single power law (i.e. F_opt < F_M,opt). The minimum expected optical flux corresponds to the case in which ν_c ∼ ν_X or ν_c > ν_X, and in this case we define the minimum allowed flux

$$\begin{array}{c}\hfill F_{m,\mathrm{opt}}=F_X{({\nu }_{\mathrm{opt}}/{\nu }_X)}^{-(\beta -1/2)}.\end{array}$$

Therefore, whenever F_m,opt < F_opt < F_M,opt, the optical and X-ray emission are compatible with a synchrotron spectrum with ν_opt < ν_c < ν_X. In the case where ν_opt < ν_c < ν_X, we can write

$$\begin{array}{c}\hfill F_{\mathrm{opt}}=F_X{\left(\frac{{\nu }_c}{{\nu }_X}\right)}^{-\beta }{\left(\frac{{\nu }_{\mathrm{opt}}}{{\nu }_c}\right)}^{-(\beta -1/2)}.\end{array}$$

Two other possibilities are:

1. F_opt > F_M,opt, and

2. F_opt < F_m,opt.

In the first case, there is no way to justify the optical and X-ray emission with a single synchrotron spectrum. In contrast, the second case can be explained if we assume that ν_X < ν_c and ν_opt < ν_m < ν_X. We recall that F_ν ∝ ν^1/3 for ν < ν_m, but none of the analysed GRBs has an optical spectral index β_opt < 0. However, we had to take into account that the spectral breaks in synchrotron are not sharp, namely, the spectral slope dF_ν/dν has a smooth transition across the break. Therefore, if ν_m ∼ ν_opt and no other break is present between ν_opt and ν_X, it is still possible to measure β_opt > 0 and F_opt < F_m,opt. Hence, we conservatively considered the second case as being compatible with a single synchrotron spectrum.

In the case of F_m,opt < F_opt < F_M,opt, since F_opt, F_X and β are known, we could estimate ν_c as follows:

$$\begin{array}{c}\hfill {\nu }_c={\left(\frac{F_X}{F_{\mathrm{opt}}}\right)}^2{\left(\frac{{\nu }_{\mathrm{opt}}}{\mathrm{keV}}\right)}^{1-p}\mathrm{keV}.\end{array}$$(2)

The uncertainty of ν_c is derived through error propagation using the same formula. Based on the comparison of optical and X-ray data, we defined two sub-samples. Sample 1 comprises all the temporal bins where the optical flux satisfies the condition F_m,opt < F_opt < F_M,opt or F_opt < F_m,opt. In this case, the optical and X-ray emissions are compatible with a single synchrotron spectrum. From our sample, 19 GRBs satisfied this condition. Sample 2 comprises data where for at least one temporal bin F_opt > F_M,opt, which indicates an incompatibility with single component synchrotron origin. From our sample, 11 GRBs fell within this case. In the case of Sample 1, we derived the temporal evolution of the break frequency ν_b, which in this case corresponds to ν_c, during and after the X-ray plateau. If the number of temporal bins where both the X-ray and optical data are available was larger than four, we fitted the temporal evolution of ν_b with the power law ν_b ∝ t^−s. The adopted fitting formula is:

$$\begin{array}{c}\hfill {log}_{10}(v_b/\mathrm{keV})=A-s{log}_{10}(\mathrm{t}/\sec ).\end{array}$$(3)

The results of the fit are reported in Table 1. In our sample, there are five GRBs where s is consistent with zero within the errors, 11 where s > 0, and three where s < 0.

We can repeat the same analysis about the temporal evolution of the spectral break relaxing the assumption that the emitting particles are in slow cooling regime. All the considerations reported above still hold, but in fast cooling we obtained:

$$\begin{array}{c}\hfill F_{m,\mathrm{opt}}=F_X{({\nu }_{\mathrm{opt}}/{\nu }_X)}^{-1/2},\end{array}$$

and the break frequency would be ν_m, which can be analogously derived and corresponds to

$$\begin{array}{c}\hfill {\nu }_m={\left(\frac{F_X}{F_{\mathrm{opt}}}\right)}^{2/(p-1)}{\left(\frac{{\nu }_{\mathrm{opt}}}{\mathrm{keV}}\right)}^{1/(1-p)}\mathrm{keV}.\end{array}$$(4)

The expression of F_M,opt is unchanged since the spectral slope above the break is the same in both the slow and fast cooling regime. This means that the selection condition F_opt ≶ F_M,opt would produce the same classification of Sample 1 and Sample 2 as in slow cooling. When comparing Eqs. (2) and (4), we noticed that the break frequency derived in the slow cooling regime (ν_c, SC) and the one derived for the fast cooling regime (ν_m, FC) are connected by the following relation:

$$\begin{array}{c}\hfill {\nu }_{m,FC}\propto {\left({\nu }_{c,SC}\right)}^{\frac{1}{p-1}}.\end{array}$$

Therefore, having fitted the temporal behaviour of ν_c, SC as ∝t^−s, in the assumption of a fast cooling regime, we would obtain

$$\begin{array}{c}\hfill {\nu }_{m,FC}\propto t^{\frac{s}{1-p}}.\end{array}$$

This implies that for typical values of p > 2, if ν_c, SC decreases in time, ν_m, FC would decrease as well.

From the values of the spectral index measured in X-rays, we can derive the value of p in the two regimes ν_c < ν_X or ν_c > ν_X. For the spectral index we considered for each GRB the time-averaged value reported in Table A.3. For Sample 1, the average X-ray spectral index is 0.94 ± 0.14, corresponding to an average value of p = 1.88 ± 0.28 (p = 2.88 ± 0.56) for ν_c < ν_X (ν_c > ν_X). For Sample 2, the average X-ray spectral index is 0.75 ± 0.14, corresponding to an average value of p = 1.51 ± 0.28 (p = 2.02 ± 0.28) for ν_c < ν_X (ν_c > ν_X).

3.3. Temporal fit of X-ray and optical light curves

In order to check if the X-ray plateau has a corresponding plateau phase in the optical, we fitted both the X-ray and optical light curves of GRBs in Sample 1 and Sample 2 with an empirical broken power law in the form:

$$\begin{array}{c}\hfill F(t)=\frac{N}{{(t/t_b)}^a+{(t/t_b)}^b}.\end{array}$$

This functional form is valid if the temporal decay before the break is shallower than the one after the break. This condition is almost always satisfied for all the X-ray and optical light curves, with the exception being the optical light curve of GRB 140419A. In this case, we adopted the following functional form:

$$\begin{array}{c}\hfill F(t)=N({(t/t_b)}^a+{(t/t_b)}^b).\end{array}$$

Regarding the optical light curves, we re-binned the data points as follows: If t_i and t_f are the initial and final times of the optical data, we re-binned using a grid of 50 bins spaced logarithmically in the interval [t_i − t_f]. The fit was performed using the python function ‘curvefit’, which is based on least squares minimisation. If the fit does not converge (i.e. if one of the parameters of the model cannot be constrained), we fitted with a single power law. The results are reported in Tables A.3 and A.4 for the X-ray and optical light curves, respectively. The best-fit curves are reported in Figs. 2 and 3. In all the cases, the fit of the X-ray light curves gave a reasonably good value of χ². However, for some optical light curves, χ² ≫ 1, indicating the presence of more complex temporal structures that cannot be approximated by a simple broken power law. The cases showing the worst agreement are GRB 100621A (Greiner et al. 2013), GRB 110213A (Cucchiara et al. 2011), and GRB 100814A (Nardini et al. 2014), which clearly exhibit the superposition of a bump over the power law decay. Such behaviour is expected, for instance, in the case of a structured jet viewed slightly off axis. Indeed, as shown by Beniamini et al. (2020b), for particular combinations of structure profile and viewing angle, the X-ray light curve shows a plateau phase, while at optical wavelengths, a bump can appear in the light curve.

Fig. 2. Summary plots of the simultaneous optical-X-ray spectral analysis of Sample 1. Here the optical and X-ray fluxes are compatible with a single synchrotron spectrum. The top panel shows the X-ray (blue points) and the optical (orange points) light curves. The blue and the orange solid lines are the best-fit curves that interpolate the X-ray and optical light curves, respectively. The optical flux density extrapolated from the X-ray flux when assuming a spectral break νb = νX is shown by fm, opt, while fM, opt is the optical flux density extrapolated from the X-ray flux when assuming no spectral break. The grey vertical band indicates when fopt < fm, opt. In the middle panel, we report the evolution of the break frequency νb as a function of time (red error bars), while the blue line indicates the best fit with a power law (nb ∝ ts). In the bottom panel, we show the X-ray spectral index for each temporal bin with black error bars, while the dot-dashed line is the time-averaged value. In the same sub-plot, we also report the optical spectral index with an orange line (the yellow band is the uncertainty interval). (a) 050319, (b) 050416A, (c) 050824, (d) 051109A.

Fig. 3. Summary plots of the GRBs from Sample 2. In this case, for at least one temporal bin, the optical flux is above fM, opt, indicating that optical and X-ray data are not compatible with a single synchrotron spectrum. The symbols are the same as in Fig. 2, and we indicate the temporal bin where fopt > fM, opt with a vertical light blue band. In the bottom panel, we show the X-ray spectral index for each temporal bin with black error bars, while the dot-dashed line is the time-averaged value. In the same sub-plot, we also report the optical spectral index with an orange line (the yellow band is the uncertainty interval). (a) 050730, (b) 060614, (c) 080310, (d) 081007.

As shown in Table A.4, eight GRBs have an optical light curve compatible with a single power law, four of which with a ∈ [ − 0.8, 0.8], thus showing the presence of a plateau phase. Nearly all the other GRBs have an optical light curve compatible with a broken power law where the flatter segment has a shallow temporal slope in the interval a ∈ [ − 0.8, 0.8]. The only exceptions are GRB 091018 (with a > 0.8), GRB 081029, and GRB 150910A (with a < −0.8). Therefore 23/30 GRBs show a plateau both in the X-ray and optical bands.

For all the GRBs that show a plateau in the optical and/or in X-ray bands, we compute the X-ray-optical luminosity during the plateau as:

$$\begin{array}{c}\hfill L_{X,\mathrm{opt}}=⟨F_{X,\mathrm{opt}}⟩\times 4\pi D_L{(z)}^2,\end{array}$$

where ⟨F_X,opt⟩ is the average flux in the specific band. The duration of the plateau was approximated with the break time. Such an approximation is valid in the limit in which the initial time of the plateau satisfies the condition t_i ≪ t_b, which is usually the case.

The relation between X-ray plateau luminosity and duration is shown in Fig. 4. The points, even if quite scattered, seem to follow the Dainotti relation (Dainotti et al. 2010), which is reported in the plot with a light blue band. The width of the band was computed by taking into account the uncertainty on both the slope and normalisation of the Dainotti relation. We fitted the L_X − t_b relation with the function:

$$\begin{array}{c}\hfill {log}_{10}\frac{L_x}{{10}^{47}\mathrm{erg}{\mathrm{s}}^{-1}}=N+\lambda {log}_{10}\frac{t_b}{{10}^4\mathrm{s}}.\end{array}$$(5)

Fig. 4. Relation between the X-ray luminosity and duration of the plateau. We indicate Sample 1 and Sample 2 with black and blue points, respectively. Arrows show upper limits. Time is reported in the rest frame of the source. The black (red) dashed line is the best-fit line for Sample 1 (Sample 2) where the fitting function is log10(LX) = N + λlog10(tb). The orange dot-dashed line is the corresponding best fit for the union of Sample 1 and Sample 2. The light blue band indicates the 1σ region defined by the best-fit values found by Dainotti et al. (2010).

We show in Fig. 5 the analogous relation between plateau luminosity and duration in the optical band. Also in this case, we fitted the same power law relation of Eq. (5). The best-fit parameters of the L_X − t_p and L_opt − t_p relations are reported in Table 2. The value of the slope of the L_X − t_b relation found by Dainotti et al. (2013) is λ = −1.32 ± 0.28, which is steeper but still compatible with our values. For the L_opt − t_b relation, Dainotti et al. (2020) report λ = −1.12 ± 0.26, which is significantly steeper than our values. We note, however, that in Dainotti et al. (2020), L_X is defined slightly differently as the luminosity at the end of the plateau.

Fig. 5. Relation between the optical luminosity and duration of the plateau. We indicate Sample 1 and Sample 2 with black and blue points, respectively. Time is reported in the rest frame of the source. The black (red) dashed line is the best-fit line for Sample 1 (Sample 2) where the fitting function is log10(Lopt) = N + λlog10(tb). The orange dot-dashed line is the corresponding best fit for the union of Sample 1 and Sample 2.

In addition to the L_X − t_p relation, previous works have shown that the properties of the X-ray plateau could correlate with the prompt emission energetics of the burst (Margutti et al. 2013; Dainotti et al. 2017; Tang et al. 2019). The investigation of these relations could give further details about the origin of the plateau. Indeed, as pointed out by Beniamini et al. (2020a), in the case of a plateau being the result of a structured jet, the X-ray plateau luminosity should correlate with the plateau duration as $L_p\propto E_{\text{iso}}{t}_p^{-1}$, where E_iso is the isotropic equivalent energy. Although, such a relation does not necessarily hold for the energy injection scenario, where the plateau energy is tightly connected to the additional energy transferred to the external shock but not to the energy released during the prompt phase. Therefore, we tested the GRBs analysed in this work to determine whether there are significant correlations between the isotropic energy E_iso and the plateau luminosity, both in the optical and X-rays bands. The correlation between plateau luminosity in X-ray band and E_iso has a Pearson coefficient of C_P = 0.58 (p-value = 7 × 10⁻⁴), while the correlation between plateau luminosity in the optical band and E_iso has a Pearson coefficient of C_P = 0.45 (p-value = 0.06). These results show that the L_X − E_iso relation is more correlated compared to the L_opt − E_iso one. These two relations are shown in panels (a) and (b) of Fig. 6.

Fig. 6. Correlations between X-ray and optical luminosities during the plateau, isotropic energy and plateau duration. In panel a (panel b) we show the correlation between the X-ray (optical) plateau luminosity and the isotropic energy Eiso. In panel c (panel d) we show the relation between luminosity and duration of the X-ray (optical) plateau, where the luminosity is corrected by the isotropic energy.

We also tested the correlation between L_X/E_iso and t_b, X as well as L_opt/E_iso and t_b, opt. The L_X − t_b, X relation has a Pearson correlation coefficient of C_P = −0.71 (p-value = 6.6 × 10⁻⁶), while the L_X/E_iso − t_b, X relation has C_P = −0.64 (p-value = 9.5 × 10⁻⁵). In the optical band, the L_opt − t_b, opt relation has a Pearson correlation coefficient of C_P = −0.25 (p-value = 0.32), while the L_opt/E_iso − t_b, opt relation has C_P = −0.41 (p-value = 0.08). Also in this case, the L_X/E_iso − t_b, X relation is much more correlated than the L_opt/E_iso − t_b, opt one, as evident from panels (c) and (d) of Fig. 6.

Finally, for the sub-sample of GRBs with a plateau in both the optical and X-ray bands, we show in Fig. 7 the relation between the duration of the plateau in the optical and X-rays bands. There is a clear indication that optical plateaus tend to be longer than X-ray plateaus, which is in line with our findings on ν_opt being smaller than ν_c in most of the cases. Among all the analysed GRBs with both an X-ray and optical plateau, only a few have t_b, X ≃ t_b, opt, implying that the majority show a chromatic plateau. There are two particular cases, GRBs 110715A and 080413B, which are in the bottom-right corner of the plot, indicating an optical break of the plateau at much later times with respect to the X-ray break. This behaviour seems to significantly deviate from the average trend. But based on Figs. A.8a and A.4a, it is possible that a simple broken power law does not sufficiently describe the temporal structure of the optical light curve, and therefore, a possibly biased estimation of the optical break could be given. Apart from these two exceptions, t_b, X and t_b, opt show a correlation. Excluding GRBs 110715A and 080413B, we derived a Pearson correlation coefficient of C_P = 0.78 and a two-tailed p-value of 3.4 × 10⁻⁴.

Fig. 7. Relation between the duration of the X-ray and optical plateau. We indicate Sample 1 and Sample 2 with black and blue points, respectively. Time is reported in the rest frame of the source. The blue dashed line corresponds to tb,opt = tb,X.

4. Discussion

The results of our analysis show that in 19 out of 30 GRBs, the plateau has an optical-to-X-ray spectrum fully consistent with synchrotron emission from a single population of shock-accelerated electrons (Sample 1). The comparison of the temporal properties of the X-ray and optical plateaus further confirms this interpretation. Indeed, the fact that the optical flux densities lie within the allowed range of values extrapolated from X-ray fluxes when assuming a single synchrotron spectrum allows us to infer ν_opt < ν_c < ν_X. This condition implies a slower evolution of the optical plateau than the X-ray one, which is in agreement with our findings that t_b, o > t_b, X in most of the GRBs belonging to Sample 1 (see Fig. 7). In order to test whether the transition of ν_c across the optical band occurs at the same time as a steepening of the optical light curve, we estimated the time t^* for which ν_c(t^*) = ν_opt and compared it with t_b, opt, where the values of the last one are reported in Table A.4. The value and the uncertainty of t^* was estimated by inverting the Eq. (3), namely:

$$\begin{array}{c}\hfill t^{\ast }={\left(\frac{{10}^N}{{\nu }_{\mathrm{opt}}/\mathrm{keV}}\right)}^{1/s}\mathrm{s}.\end{array}$$

In Table 3, we compare the t^* and t_b, opt of all the GRBs of Sample 1 that have s > 0. The only case where t^* and t_b, opt are not compatible is GRB 090618. In all other cases, we have three possibilities: (1) We have an estimate of both t^* and t_b, opt, and they are compatible. (2) Only a lower limit on t^* is available but still compatible with t_b, opt. (3) No break in the optical is observed, possibly indicating that the break may occur later than the period when the data was collected. The second and third cases are inconclusive, and only for the first case (GRBs 061121, 080413B, and 111228A) we were able to conclude that t^* ≃ t_b, opt, indicating that the optical temporal break is in agreement with the transition of the cooling frequency across the optical band.

If the condition ν_opt < ν_c < ν_X is satisfied, no spectral evolution is expected in the X-rays. We tested the presence of spectral evolution by fitting the X-ray photon index temporal trend with a constant. In Table A.3, we report the average photon index (the average was taken over the whole available data) and the p-value of the fit. Of the 19 cases in Sample 1, six GRBs have a p-value less than 0.05, showing significant deviation from the constant trend. For three of these cases (GRBs 060729, 061121, and 090618), the deviation is given by the hardening of the X-ray spectrum at late times. For remaining three cases, GRBs 100621A, 110213A, and 140419A, there is evidence of spectral softening. For five of the six cases, the exception being GRB 140419A, Table 1 shows that the derived ν_c decreases in time. In this scenario, a softening of the X-ray spectrum is expected, as F_ν ∝ ν^{−(p − 1)/2} for ν < ν_c and F_ν ∝ ν^−p/2 for ν > ν_c. Since the synchrotron spectral shape has smooth transitions between two adjacent power law segments, a temporal decrease of ν_c would produce gradual spectral softening. However, a trend of a decreasing ν_c is not in agreement with the evidence of spectral hardening. In this regard, we point out that in the derivation of the temporal evolution of ν_c, we assumed the spectral slopes above and below ν_c to be constant in time. If there is an additional process that induces an intrinsic variation of the spectral slope, the estimation of the temporal trend of ν_c can be biased. One possible cause for the spectral hardening observed in the X-rays could be the additional cooling of electrons by SSC emission. Indeed, as shown by Beniamini et al. (2015), for instance, the emergence of the SSC component can cause a transition from the spectral slope p/2 (above the cooling frequency) to a harder value. An analogous effect can be obtained if the synchrotron emission is still dominant but the Klein-Nishina (KN) effect is non-negligible (e.g. Nakar et al. 2009). For typical micro-physical parameters, Nakar et al. (2009) estimated that when including the KN effect, not only is the temporal evolution of the cooling frequency modified (${\nu }_c\propto t^{-\frac{8-3p}{4-p}}$), but the spectral index in the X-rays (for ν_X > ν_c) would also change, from p/2 to 0.75(p − 1), thereby causing a slight hardening in the X-rays. Hence, we conclude that the results shown for the evolution of the cooling frequency are fully reliable for the GRBs that do not show significant spectral evolution in the X-rays. Finally, amongst the 11 cases of Sample 2, six have a p-value < 0.05 (GRBs 050730, 080310, 100814A, 110715A, 120404A, and 150910A), but none of them show a clear trend that points towards a softening or a hardening of the X-ray spectrum.

4.1. Modifications to the standard forward shock scenario

The consistency of the optical and X-ray data with a single spectrum can be interpreted as an indication that both X-ray and optical photons originate from the same emission region. Moreover, the process responsible for the X-ray plateau should also explain the observed evolution of the cooling frequency. In the standard scenario of a forward shock decelerating through the circumburst medium, the predicted temporal evolution of the characteristic synchrotron frequencies is (Granot & Sari 2002):

$$\begin{array}{c}\hfill {\nu }_c\propto {ϵ}_{\mathrm{B}}^{-3/2}{n}_0^{-1}{E}_{52}^{-1/2}{t}_{\text{days}}^{-1/2}\end{array}$$(6)

for an interstellar medium (ISM) with uniform particle density, n₀, and

$$\begin{array}{c}\hfill {\nu }_c\propto {ϵ}_{\mathrm{B}}^{-3/2}{E}_{52}^{1/2}{t}_{\mathrm{days}}^{1/2}\end{array}$$(7)

for a medium with a wind-like density profile. With E₅₂, we indicate the isotropic energy in units of 10⁵² erg. The flux density is expected to decline as F_ν ∝ t^−a, and the predicted values are in a slow cooling regime (e.g. see Zhang et al. 2006):

$$\begin{array}{c}\hfill a=\frac{3}{4}(p-1)\text{for}\nu <{\nu }_c,a=\frac{3p-2}{4}\text{for}\nu >{\nu }_c\end{array}$$

for the ISM scenario and

$$\begin{array}{c}\hfill a=\frac{3p-1}{4}\text{for}\nu <{\nu }_c,a=\frac{3p-2}{4}\text{for}\nu >{\nu }_c\end{array}$$

for the wind scenario. Since the values of p are likely above two, the standard picture predicts afterglow light curves no flatter than t^−3/4. Therefore, even if the observed temporal evolution of ν_c can be compatible with the temporal behaviour specified in Eqs. (6) and (7), the observed temporal slope of the X-ray and optical plateaus are incompatible with the standard forward shock scenario.

Modifications to the standard picture have been proposed in the literature (Granot et al. 2006; Ioka et al. 2006; Panaitescu et al. 2006; Misra et al. 2021) and invoke the possible temporal evolution of the shock micro-physical parameters, such as ϵ_B. In the specific, assuming a temporal evolution in the form ϵ_B ∝ t^μ and a circumburst medium with a density profile ρ ∝ r^−k, the cooling frequency evolves as ν_c ∝ t^−s, with

$$\begin{array}{c}\hfill s=\frac{4+12\mu -3k(1+\mu )}{2(4-k)}.\end{array}$$

Moreover, when considering the ordering ν_opt < ν_c < ν_X, the optical and X-ray flux densities are expected to evolve as F_opt ∝ t^a_o and F_X ∝ t^a_X, with

$$\begin{array}{c}\hfill a_o=\frac{1}{2}(\mu +\frac{k}{k-4})+\frac{p-1}{2}\left(\frac{4+12\mu -3k(1+\mu )}{2(k-4)}\right)\end{array}$$

and

$$\begin{array}{c}\hfill a_X=\frac{1}{2}(\mu +\frac{k}{k-4})+\frac{p}{2}\left(\frac{4+12\mu -3k(1+\mu )}{2(k-4)}\right).\end{array}$$

In order to check the validity of this scenario, we compared the expected values of a_o, a_X, and s with the observed values. In particular, we searched for a combination of μ and k such that the following relations are simultaneously satisfied:

$$\begin{array}{c}\hfill a_{o,th}\in [a_{o,\mathrm{obs}}-\mathrm{\Delta }{a}_o,a_{o,\mathrm{obs}}+\mathrm{\Delta }{a}_o],\end{array}$$

$$\begin{array}{c}\hfill a_{X,th}\in [a_{X,\mathrm{obs}}-\mathrm{\Delta }{a}_X,a_{X,\mathrm{obs}}+\mathrm{\Delta }{a}_X],\end{array}$$

$$\begin{array}{c}\hfill s_{\mathit{th}}\in [s_{\mathrm{obs}}-\mathrm{\Delta }s,s_{\mathrm{obs}}+\mathrm{\Delta }s],\end{array}$$

where we indicate with X_th, X_obs, and ΔX the expected value, the observed value, and the error of the quantity X, respectively. The intersection of these three conditions defines a region in the k − μ plane that is in agreement with observations. We verified for all the GRBs of Sample 1 that a combination of k and μ that satisfies the above conditions does not exist.

4.2. The energy injection scenario

Another solution invokes the presence of additional energy injected into the forward shock at late times. A single emitting region is compatible with the energy injection scenario, where additional energy is transferred to the external shock due to the late-time activity of the central engine. In this scenario, the dynamical evolution of the blast wave is determined by the rate of energy injection and the efficiency of conversion from injected energy to jet kinetic energy. Hence, the deceleration in the ISM is less severe and the flux drop is shallower. The injection of energy modifies the dynamical evolution of the external shock, but it does not have an influence on the dominating radiative mechanisms responsible for the dissipation of the particles’ energy. Therefore, the spectral energy distribution should be the same of that in the standard scenario of particles dissipating through synchrotron radiation. We note that if the energy injection has an impact only on the blast wave dynamics, the plateau should be achromatic (e.g. Fan & Piran 2006). This would imply that the temporal break corresponding to the transition from plateau to post-plateau phase should be the same in the optical and X-rays bands. Although, as pointed out before, this is not the case if a spectral break is present between the two bands.

When adopting the standard prescription of a injected luminosity in the form

$$\begin{array}{c}\hfill L_{\mathrm{inj}}\propto t^{-q},\end{array}$$

the temporal slope of the afterglow light curve at any frequency ν_obs can be predicted depending on the position of ν_obs relative to the synchrotron characteristic frequencies (ν_m and ν_c). Moreover the temporal behaviour of these frequencies are modified consequently and read as

$$\begin{array}{c}\hfill {\nu }_m\propto t^{-(2+q)/2},{\nu }_c\propto t^{(q-2)/2}\end{array}$$(8)

for the ISM scenario (jet propagating into a medium with constant density) and

$$\begin{array}{c}\hfill {\nu }_m\propto t^{-(2+q)/2},{\nu }_c\propto t^{(2-q)/2}\end{array}$$(9)

for the wind scenario (jet propagating into a medium with density n ∝ r⁻²). Therefore, if q is less than two (q is greater than two) the energy injection model predicts a decreasing cooling frequency in the ISM scenario (wind scenario).

One of the possible sources of energy responsible for the jet refreshing is the long-lived, highly magnetised neutron star (magnetar) left after the production of the GRB (Dai & Lu 1998; Zhang & Mészáros 2001; Dall’Osso et al. 2011). The magnetar loses rotational energy through spin-down radiation. The associated released luminosity in the standard scenario of a rotating magnetic dipole depends on the rotational frequency as L_sd ∝ Ω⁴. The temporal evolution of the spin-down luminosity is:

$$\begin{array}{c}\hfill L_{\mathrm{sd}}(t)=\frac{L_0}{{[1+\frac{t}{\tau }]}^2}.\end{array}$$

This relation can be further extended by including deviations from the standard picture of simple dipole radiation. This modification leads to a spin-down luminosity in the form L_sd ∝ Ω^{4 − 2α}, where α is related to the braking index n as n = 3 − 2α and 0 < α < 1. The corresponding temporal behaviour of the spin-down luminosity is:

$$\begin{array}{c}\hfill L_{\mathrm{sd}}(t)=\frac{L_0}{{[1+(1-\alpha )\frac{t}{\tau }]}^{\frac{2-\alpha }{1-\alpha }}}.\end{array}$$

For t ≫ τ, the spin-down luminosity evolves as L_sd ∝ t^−q, with

$$\begin{array}{c}\hfill q=\frac{2-\alpha }{1-\alpha }⩾2,\end{array}$$(10)

and the forward shock afterglow emission dominates. Indeed, following Dall’Osso et al. (2011), the luminosity evolution in the relativistic external shock can be obtained from the balance between radiative losses and energy injection from the spinning-down magnetar. Fitting this model on the observed X-ray light curves provided a very good description of the plateau and the post-plateau phases, giving reasonable values of the magnetic field strength and spin period (e.g. Dall’Osso et al. 2011; Bernardini et al. 2012; Stratta et al. 2018). In this scenario, the post-plateau afterglow spectral properties follow the standard forward shock prescriptions and is in good agreement with the majority of the GRBs in Sample 1. Since the majority of the GRBs in Sample 1 (16/19) show a preference for ν_c decreasing in time after the plateau phase, even in the magnetar scenario we can apply Eqs. (6) and (7) and find agreement with a forward shock propagating in an ISM medium. For the three remaining GRBs, ν_c increases with time with a slope compatible with t^0.5, indicating a preference for a wind-like density profile.

If the X-ray plateau is a direct consequence of additional energy being injected into the forward shock, it would also impact the ratio between the jet kinetic energy E_kin (inferred by the afterglow light curve) and the prompt isotropic energy E_iso. The ratio ϵ_pr = E_iso/E_kin represents the prompt emission efficiency (Beniamini et al. 2016), and if the injection of energy in the forward shock tends to increase E_kin, then the inferred ϵ_pr should be lower with respect to average cases where no energy injection operates. In order to verify this point, we followed the method adopted by D’Avanzo et al. (2012), which uses the X-ray luminosity at 24 hr after the burst (L_X, 24 h) as a direct proxy of E_kin. In Fig. 8, we show the relation between L_X, 24 h and E_iso, comparing our GRB sample with the average relation reported in Fig. 2 by D’Avanzo et al. (2012), where the authors consider a complete sample of GRBs detected by Swift with no specific selection regarding the presence of an X-ray plateau. Our plot shows that our sample of GRBs is characterised by a ratio L_X, 24 h/E_iso larger than the average, which can be an indication of a larger ratio E_kin/E_iso and therefore favour the energy injection scenario. We note that the latter implication is based on the adoption of a specific set of micro-physical parameters of the external shock. Therefore, the conclusion that the ratio E_kin/E_iso is larger for our sample holds in the assumption that all GRBs (with and without a plateau) share similar micro-physical parameters.

Fig. 8. Comparison between the isotropic equivalent energy Eiso and the X-ray luminosity at 24 h after the trigger. The blue solid line and the yellow band represent the best fit and 3σ dispersion, respectively, of the same relation reported by D’Avanzo et al. (2012) for a complete sample of GRBs observed by Swift.

4.3. The plateau as a result of the jet structure

The HLE from a structured jet also provides a viable explanation for the presence of a plateau phase in the X-ray light curves of GRBs. In the approach of Oganesyan et al. (2020), the photons emitted during the prompt emission at large angles with respect to the jet core arrive to the observer at late times and are less Doppler boosted. The duration and temporal slope of the plateau phase depend on the jet structure, radius, and the Lorentz factor of the dissipation site. In the approximation of an instantaneous prompt emission and in the limit of an on-axis observer, there is a bi-univocal relation between the arrival time t_obs of the photon and θ, where θ is the polar angle between the jet axis and the patch of the shock front from which the photon departed. Namely, the flux observed at time t_obs corresponds to the contribution from a jet ring at a polar angle θ(t_obs). This implies that, if the energy spectrum in the shock comoving frame is the same along the jet ring, then the observed spectrum at each time should reflect the same spectral shape of the prompt emission. Though, since the Doppler boosting decreases at higher latitudes, the observed spectral peak should decrease with time as well. Therefore, also in the HLE scenario, we expect an optical and X-ray emission compatible with a single synchrotron component whose characteristic frequencies (ν_m and ν_c) decrease in time due to the Doppler effect. The specific rate of the temporal decrease of the characteristic frequencies depends on the jet structure. We note here that HLE would be dominant over the afterglow from forward shock emission throughout the light curve, as it would still be compatible with a single synchrotron component.

Another approach based on the HLE from a structured jet is the one adopted by Beniamini et al. (2020a,b), where the photons responsible for the plateau emission are produced in the deceleration phase of the forward shock. In particular, in Beniamini et al. (2022), the authors show the predictions of the temporal evolution of the cooling frequency. As shown in Fig. 1 and in Table 1 of Beniamini et al. (2022), the model predicts

$$\begin{array}{c}\hfill -1<s=-\frac{3k-4}{2}<2,\text{for}0<k<2\end{array}$$(11)

for the pre-deceleration and the angular structure dominated phase, while it predicts

$$\begin{array}{c}\hfill -1/2<s=-\frac{3k-4}{8-2k}<1/2,\text{for}0<k<2\end{array}$$(12)

for the post deceleration phase. As before, ρ ∝ r^−k. These ranges of s are fully consistent with the ones reported in Table 1. In this class of models, the plateau appears when the jet is observed slightly off axis, namely, at θ_obs ≳ θ_j, where θ_j is the opening angle of the jet core. In Fig. 9, the plateau sample shows no significant deviation from the average E_iso − E_p correlation, which is largely consistent with the condition θ_obs ∼ θ_j.

Fig. 9. Ep − Eiso plane. We report in orange the GRB sample selected in our work, see Table A.1. The blue points were taken from Nava et al. (2012), and the grey region is the corresponding 3σ band.

4.4. Alternative scenarios

In contrast to the energy injection scenario, the emission from the forward shock can still produce a plateau phase if the blast wave is stratified in velocity (Rees & Mészáros 1998). More specifically, if the faster component of the ejecta is launched later, it catches up to the slower one, resulting in a refreshed shock wave. If the stratification is given by a smooth distribution in the Lorentz factor, the effect is a flattening of the light curve, while in the limit of two distinct shells launched at different speeds, the light curve could show a re-brightening followed again by standard forward shock decay. However, we notice that in this scenario, the broad-band emission should also be compatible with a single synchrotron component.

Provided that the forward shock radiative efficiency is sufficiently low, the reverse shock emission can dominate during the plateau phase as well. Uhm & Beloborodov (2007) demonstrated that the reverse shock emission from a stratified blast wave can easily produce a flat segment in the X-ray light curve, without assuming any injection of energy in the forward shock (see also Genet et al. 2007). Moreover the reverse shock emission is sensitive to the change of density in the external medium, which corresponds to a re-brightening of the optical light curve, as observed for several GRBs in Sample 1 (GRBs 100621A, 110213A) and in Sample 2 (GRB 100814A).

Another interesting scenario considers the possibility that the plateau phase occurs in the coasting phase of jet propagation, namely, before blast wave deceleration (Granot et al. 2006; Shen & Matzner 2012). Following this idea, Dereli-Bégué et al. (2022) claim that the plateau can still be explained in the classical fireball model, provided that the jet is propagating in a wind environment and with a rather low bulk Lorentz factor value (approximately a few tens). This last condition ensures the jet to have a long enough coasting phase and a corresponding light curve whose temporal slope is close to the ones observed during the plateau. However, in order for the jet dissipation radius to be greater than the photospheric radius, such low values of the Lorentz factor would require an isotropic luminosity that is not too large and a minimum variability time scale of few seconds. Moreover, the transparency requirement is satisfied if rather high values of ϵ_B (> 0.1) are assumed. Therefore, this scenario also requires further investigation.

4.5. Remarks about Sample 2

For the GRBs belonging to Sample 2, optical and X-ray data are incompatible with a single synchrotron spectrum in at least one temporal bin. In two events (GRB 100814A and GRB 150910A), the optical light curve has a shape that differs substantially from that of the X-ray light curve. In a few cases, though, the two light curves behave in a similar way, yet the optical emission appears to be a few times above the maximum extrapolated from the X-ray flux (e.g. GRB 060614, GRB 081007, and GRB 100418A). Amongst the GRBs in Sample 1, the case of GRB 100621A shows several peculiarities (Greiner et al. 2013). The optical light curve clearly shows a bump around 5 × 10³ − 2 × 10⁴ s, while the X-ray light curve is well fitted by a broken power law. Moreover, the X-ray spectrum shows an evident softening at late times. Such a softening imposes the introduction of an additional spectral break between the optical and X-rays bands in order to explain the broad-band emission with a single synchrotron component. If this spectral break is identified with ν_m, the observations would imply a ν_m that increases with time. This increase is hardly explained in the framework of standard forward shock theory, even when including an energy injection term (see Eqs. (8) and (9)). Hence, even if GRB 100621A satisfies the conditions to be in Sample 1, the shape of the optical light curve and the X-ray spectral softening are hardly explained in the context of standard synchrotron emission from a single emission zone.

When the optical light curve shows substantial differences with the X-ray one, a single synchrotron-emitting region cannot explain the observed broad-band spectrum. One viable explanation is that the HLE associated with the prompt phase and the external shock emission from a structured jet can simultaneously contribute to the broad-band spectrum. The HLE can dominate in the X-rays, while the optical band would be dominated by external shock. If this last option is valid, it means that the two components cannot have the same spectral energy distribution. Otherwise, one of them would dominate in both of the bands, which is not the case. Moreover, if the component that dominates in the optical band is characterised by a non-thermal spectrum, we should expect an optical spectral index softer than the X-ray one. Indeed, in order to have an excess flux in the optical band, the dominating spectral component in the band should have a steeper spectral slope. If in the optical band F_ν ∼ ν^−β_opt and in the X-ray band F_ν ∼ ν^−β_X, this would translate to β_opt > β_X. Nonetheless, when comparing the optical spectral slopes with the X-ray spectral index for Sample 2, we verified that they have compatible values. A possible solution would thus be that the optical flux is dominated by a component with a non-thermal spectrum and a cut-off at energies ν_opt < ν_{cut − off} < ν_X or, alternatively, by a thermal component with a characteristic temperature kT ∼ hν_opt.

5. Conclusions

The origin of the plateau phase in the X-ray light curve of GRBs is still a matter of debate. In this work, we analysed a complete sample of GRBs with simultaneous optical and X-ray data during the plateau phase in order to shed light on its physical origin. We performed a time-resolved spectral analysis in the X-rays, and we compared X-ray and optical data to verify if they are compatible with a single synchrotron spectrum. While the majority of the cases show compatibility with a single component synchrotron origin of the multi-band emission, we collected evidence that some GRBs (approximately one-third of the entire sample) are incompatible with the standard forward shock emission from a single dissipation zone.

For the sample of GRBs compatible with a single-zone emission, we derived the temporal evolution of the cooling frequency, and we compared it with the predictions from several models. We show that the majority of the GRBs show a cooling frequency that decreases in time. For the GRBs that show a plateau both in the X-ray and in the optical bands (19 over 30), we found that the optical plateau has typically a longer duration than the X-ray one.

Moreover, we verified that the transition of the cooling frequency across the optical band is compatible with a simultaneous steepening of the optical light curve. We find that both an energy injection model and the scenario of HLE from a structured jet predict a temporal decay slope in the optical-X-ray and a trend of the cooling frequency compatible with the observations, since both scenarios assume classic synchrotron emission. However, while this emission would come from the afterglow (external shock) in the energy injection scenario, the HLE emission model would imply that everything comes from the prompt emission region and the external shock would always remain sub-dominant. A model in which the plateau is produced by the external shock viewed off axis, as the one proposed by Beniamini et al. (2022), would be degenerate with energy injection in the framework of our study. Moreover, the results for our plateau sample show that the ratio between the X-ray luminosity at 24 h and the isotropic equivalent energy is on average larger than in other GRBs, possibly favouring the energy injection scenario.

Concerning the second sample of GRBs not compatible with a single synchrotron spectrum, the optical emission lies above the extrapolation inferred from the X-ray analysis, at least in some time intervals during the plateau. Neither the energy injection model nor the HLE model alone can account for this behaviour. Such a result necessarily requires going beyond the standard modelling that we used here. Allowing for a more complex spectral model and/or the co-existence of two (or more) emission regions is likely to provide a satisfactory solution. A confirmation of this ‘mixed’ scenario requires further investigation through the detailed modelling of the different emission components.

Acknowledgments

We acknowledge financial contribution from the agreement ASI-INAF n.2017-14-H.O. GS and DAK acknowledge the support by the State of Hesse within the Research Cluster ELEMENTS (Project ID 500/10.006). DAK acknowledges support from Spanish National Research Project RTI2018-098104-J-I00 (GRBPhot). MB and GO acknowledge financial support from the AHEAD2020 project (grant agreement n. 871158). AR acknowledges support from the INAF project Premiale Supporto Arizona & Italia. SD acknowledges funding from the European Union’s Horizon2020 research and innovation program under the Marie Skłodowska-Curie (grant agreement No.754496). Sadly, one of the authors, David Alexander Kann, passed away before the publication of this paper. His premature departure is a deep loss for the astronomical community. We spend these words to remember his outstanding achievements, his unique insight, and invaluable contributions in the fields of gamma-ray bursts and multi-messenger astronomy.

References

Appendix A: The role of the Y parameter

Following Sari & Esin (2001), here we estimate how the Y parameter adopted in eq. 1 can influence the ratio ν_c/ν_m and hence determine the cooling regime. The Y parameter is defined as the ratio between the SSC and synchrotron luminosity:

$$\begin{array}{c}\hfill Y=\frac{L_{\mathrm{SSC}}}{L_{\mathrm{syn}}}=\frac{{\eta }_e{ϵ}_e}{{ϵ}_B(1+Y)},\end{array}$$

where

$$\begin{array}{c}\hfill {\eta }_e=\{\begin{array}{cc}1,\hfill & \text{fast cooling}\hfill \\ {\left(\frac{{\gamma }_c}{{\gamma }_m}\right)}^{2-p},\hfill & \text{slow cooling}\hfill \end{array},\end{array}$$

from which

$$\begin{array}{c}\hfill Y=\{\begin{array}{cc}\frac{{\eta }_e{ϵ}_e}{{ϵ}_B},\hfill & \frac{{\eta }_e{ϵ}_e}{{ϵ}_B}\ll 1\hfill \\ {\left(\frac{{\eta }_e{ϵ}_e}{{ϵ}_B}\right)}^{1/2},\hfill & \frac{{\eta }_e{ϵ}_e}{{ϵ}_B}\gg 1\hfill \end{array}.\end{array}$$

Next, we define ν_c/ν_m = (γ_c/γ_m)² = R. Since R depends on Y, which in turn can depend on R, we check a posteriori the value of R, depending on the value of Y. As reference, we assume ϵ_e = 0.1 and ϵ_B = 10⁻³ as fiducial values.

1. If we assume that we are in slow cooling, then R > 1. We can have:

   • (a)

     Y ≪ 1, namely η_eϵ_e/ϵ_B ≪ 1, which gives $R\ll {\left(\frac{{ϵ}_e}{{ϵ}_B}\right)}^{\frac{2}{p-2}}$. This condition, for p > 2, is compatible with the value of R derived from eq. 1 (R∼10⁵); hence, it is not in contradiction with a slow cooling regime. In other words, slow cooling and negligible SSC can co-exist in the condition $1<R\ll {\left(\frac{{ϵ}_e}{{ϵ}_B}\right)}^{\frac{2}{p-2}}$.

   • (b)

     Y ≫ 1, namely η_eϵ_e/ϵ_B ≫ 1, leading to $R\gg {\left(\frac{{ϵ}_e}{{ϵ}_B}\right)}^{\frac{2}{p-2}}\gg 1$. Though, if the expression of $Y=\sqrt{R^{\frac{2-p}{2}}\frac{{ϵ}_e}{{ϵ}_B}}$ is inserted into eq. 1, we obtain

     $$\begin{array}{c}\hfill R\sim {\left({10}^5\frac{{ϵ}_B}{{ϵ}_e}\right)}^{\frac{1}{2-p/2}}\ll {\left(\frac{{ϵ}_e}{{ϵ}_B}\right)}^{\frac{2}{p-2}},\end{array}$$

     which creates a contradiction. This means that, for the assumed values of ϵ_e and ϵ_B, the condition Y ≫ 1 is incompatible with a slow cooling regime.

2. If we assume a fast cooling regime, η_e = 1 and $Y=\sqrt{{ϵ}_e/{ϵ}_B}=10$. Substituting Y in eq. 1, we get R ∼ 10³, which is in contradiction with the assumption of a fast cooling.

To summarise, with the assumed values of ϵ_e and ϵ_B, the only acceptable combination is Y ≪ 1 and a slow cooling regime. More generally, the value of R can be found as a function of ϵ_e and ϵ_B, solving numerically the implicit function defined by eq. 1, where R = R(ϵ_e, ϵ_B, Y(R, ϵ_e, ϵ_B)). The result of this numerical solution is shown in Fig. A.1, where we fixed t_d = 1 and n = 1 cm⁻³. The plot shows that, for the assumed fiducial values ϵ_e = 0.1 andϵ_B = 10⁻³, the system is in the slow cooling regime. Fig. A.2 shows the same dividing line from Fig. A.1 but at three different times: t = 10⁻² d, t = 10⁻¹ d, and t = 1 d. It is evident that at early times, the region of the [ϵ_e, ϵ_B] plane, which gives fast cooling, is larger. Though for the assumed fiducial values ϵ_e = 0.1 and ϵ_B = 10⁻³, at early times, the system is still in the slow cooling regime. We note that the results shown in this section do not rule out the fast cooling regime but rather demonstrate that slow cooling is more plausible since there exist more combinations of ϵ_e and ϵ_B that give ν_c > ν_m. For this reason, since fast cooling cannot be rejected, we interpret our data in both cooling regimes in section 3.2.

Fig. A.1. Contour plot of the value of R = νc/νm as a function of ϵe and ϵB. The dashed line corresponds to the combination of ϵe and ϵB, which gives R = 1, and it divides the plane between the slow cooling (bottom left) and the fast cooling (top right) regimes. The region on the right corresponds to the fast cooling regime.

Fig. A.2. Lines in the [ϵe, ϵB] plane where R = νc/νm = 1. Each line was computed at different times. The region of the plane to the right of each line corresponds to the combinations of [ϵe, ϵB] that give a fast cooling regime.

Fig. A.3. Sample 1 - continued.

Fig. A.4. Sample 1 - continued.

Fig. A.5. Sample 1 - continued.

Fig. A.6. Sample 1 - continued.

Fig. A.7. Sample 2 - continued.

Fig. A.8. Sample 2 - continued.

All Tables

All Figures

Fig. 1. Schematic representation of the method used to compare X-ray and optical data for each temporal bin. The minimum and maximum optical fluxes allowed by the standard afterglow model assuming that there is a single break νopt < νb < νX are indicated by Fm,opt and FM,opt, respectively. We derived FX, Fopt, and the X-ray spectral slope from spectral analysis. Therefore, if Fm,opt < Fopt < FM,opt, we then considered that the optical counterpart is consistent with a single synchrotron spectrum, and we could derive νb. The blue line has no break, the orange line has νb = νX, and the green line has a break between the optical and X-ray bands. We note that our assumption here is that νm < νopt, and it is supported by the observation of βopt > 0 for all the analysed events (see Sect. 3.2) that exclude the case νb = νm.

In the text

Fig. 2. Summary plots of the simultaneous optical-X-ray spectral analysis of Sample 1. Here the optical and X-ray fluxes are compatible with a single synchrotron spectrum. The top panel shows the X-ray (blue points) and the optical (orange points) light curves. The blue and the orange solid lines are the best-fit curves that interpolate the X-ray and optical light curves, respectively. The optical flux density extrapolated from the X-ray flux when assuming a spectral break νb = νX is shown by fm, opt, while fM, opt is the optical flux density extrapolated from the X-ray flux when assuming no spectral break. The grey vertical band indicates when fopt < fm, opt. In the middle panel, we report the evolution of the break frequency νb as a function of time (red error bars), while the blue line indicates the best fit with a power law (nb ∝ ts). In the bottom panel, we show the X-ray spectral index for each temporal bin with black error bars, while the dot-dashed line is the time-averaged value. In the same sub-plot, we also report the optical spectral index with an orange line (the yellow band is the uncertainty interval). (a) 050319, (b) 050416A, (c) 050824, (d) 051109A.

In the text

Fig. 3. Summary plots of the GRBs from Sample 2. In this case, for at least one temporal bin, the optical flux is above fM, opt, indicating that optical and X-ray data are not compatible with a single synchrotron spectrum. The symbols are the same as in Fig. 2, and we indicate the temporal bin where fopt > fM, opt with a vertical light blue band. In the bottom panel, we show the X-ray spectral index for each temporal bin with black error bars, while the dot-dashed line is the time-averaged value. In the same sub-plot, we also report the optical spectral index with an orange line (the yellow band is the uncertainty interval). (a) 050730, (b) 060614, (c) 080310, (d) 081007.

In the text

Fig. 4. Relation between the X-ray luminosity and duration of the plateau. We indicate Sample 1 and Sample 2 with black and blue points, respectively. Arrows show upper limits. Time is reported in the rest frame of the source. The black (red) dashed line is the best-fit line for Sample 1 (Sample 2) where the fitting function is log10(LX) = N + λlog10(tb). The orange dot-dashed line is the corresponding best fit for the union of Sample 1 and Sample 2. The light blue band indicates the 1σ region defined by the best-fit values found by Dainotti et al. (2010).

In the text

Fig. 5. Relation between the optical luminosity and duration of the plateau. We indicate Sample 1 and Sample 2 with black and blue points, respectively. Time is reported in the rest frame of the source. The black (red) dashed line is the best-fit line for Sample 1 (Sample 2) where the fitting function is log10(Lopt) = N + λlog10(tb). The orange dot-dashed line is the corresponding best fit for the union of Sample 1 and Sample 2.

In the text

Fig. 6. Correlations between X-ray and optical luminosities during the plateau, isotropic energy and plateau duration. In panel a (panel b) we show the correlation between the X-ray (optical) plateau luminosity and the isotropic energy Eiso. In panel c (panel d) we show the relation between luminosity and duration of the X-ray (optical) plateau, where the luminosity is corrected by the isotropic energy.

In the text

	Fig. 7. Relation between the duration of the X-ray and optical plateau. We indicate Sample 1 and Sample 2 with black and blue points, respectively. Time is reported in the rest frame of the source. The blue dashed line corresponds to tb,opt = tb,X.
In the text

	Fig. 8. Comparison between the isotropic equivalent energy Eiso and the X-ray luminosity at 24 h after the trigger. The blue solid line and the yellow band represent the best fit and 3σ dispersion, respectively, of the same relation reported by D’Avanzo et al. (2012) for a complete sample of GRBs observed by Swift.
In the text

	Fig. 9. Ep − Eiso plane. We report in orange the GRB sample selected in our work, see Table A.1. The blue points were taken from Nava et al. (2012), and the grey region is the corresponding 3σ band.
In the text

	Fig. A.1. Contour plot of the value of R = νc/νm as a function of ϵe and ϵB. The dashed line corresponds to the combination of ϵe and ϵB, which gives R = 1, and it divides the plane between the slow cooling (bottom left) and the fast cooling (top right) regimes. The region on the right corresponds to the fast cooling regime.
In the text

	Fig. A.2. Lines in the [ϵe, ϵB] plane where R = νc/νm = 1. Each line was computed at different times. The region of the plane to the right of each line corresponds to the combinations of [ϵe, ϵB] that give a fast cooling regime.
In the text

	Fig. A.3. Sample 1 - continued.
In the text

	Fig. A.4. Sample 1 - continued.
In the text

	Fig. A.5. Sample 1 - continued.
In the text

	Fig. A.6. Sample 1 - continued.
In the text

	Fig. A.7. Sample 2 - continued.
In the text

	Fig. A.8. Sample 2 - continued.
In the text