Unraveling Parameter Degeneracy in GRB Data Analysis

Gamma-ray burst (GRB) afterglow light curves and spectra provide information about the density of the environment, the energy of the explosion, the properties of the particle acceleration process, and the structure of the decelerating jet. Due to the large number of parameters involved, the model can present a certain degree of parameter degeneracy. In this paper, we generated synthetic photometric data points using a standard GRB afterglow model and fit them using the Markov Chain Monte Carlo (MCMC) method. This method has emerged as the preferred approach for analysing and interpreting data in astronomy. We show that, depending on the choice of priors, the parameter degeneracy can go unnoticed by the MCMC method. Furthermore, we apply the MCMC method to analyse the GRB~170817A afterglow. We find that there is a complete degeneracy between the energy of the explosion $E$, the density of the environment $n$, and the microphysical parameters describing the particle acceleration process (e.g. $\epsilon_e$ and $\epsilon_B$), which cannot be determined by the afterglow light curve alone. Our results emphasise the importance of gaining a deep understanding of the degeneracy properties which can be present in GRB afterglows models, as well as the limitations of the MCMC method. In the case of GRB 170817, we get the following values for the physical parameters: $E=8\times 10^{50}-1 \times 10^{53}$ erg, $n=7\times 10^{-5}-9\times10^{-3}$, $\epsilon_e=10^{-3}-0.3$, $\epsilon_B=10^{-10}-0.3$.


INTRODUCTION
Gamma-ray bursts (GRBs) are generated during the propagation of highly relativistic jets (see, e.g., Kumar & Zhang 2015, for a review).The gamma-ray emission results from the dissipation of kinetic or magnetic energy within the jet channel, while the subsequent multiwavelength afterglow is a consequence of the interaction between the relativistic jet and the surrounding medium (Mészáros & Rees 1997).
The afterglow phase can be generally 1 explained as synchrotron radiation from a population of electrons accelerated by the relativistic shocks associated with the decelerating jet (see, e.g.Miceli & Nava 2022, for a recent review).In the most basic model of a spherical explosion propagating through a uniform medium, the Lorentz factor and velocity of the blast wave go as Γ sh ∝ t −3/2 during the relativistic phase, and v sh ∝ t −3/2 when the blast wave becomes sub-relativistic, being t the time in the laboratory frame (Blandford & Mc-Kee 1976).Then, the afterglow emission is typically characterised by a series of power-law segments, with an observed flux density given as Fν ∝ t −α ν −β , being t and ν the time and frequency in the observer frame, and α and β the tempo-1 Very high-energy emission from GRB 190829A (H.E. S. S. Collaboration et al. 2021), GRB 190114C (Ajello et al. 2020) and GRB 180720B (Abdalla et al. 2019) are signatures of an inverse Compton component after the prompt emission.
Markov chain Monte Carlo (MCMC) methods (see e.g.Sharma 2017; Hogg & Foreman-Mackey 2018; Speagle 2019), in particular, are the most widely used strategies currently employed to model GRB observations.These techniques have been used in the GRB community to determine not only multifrequency adjustments (e.g., Dobie et al. 2018;Mooley et al. 2018b;Bright et al. 2019;Makhathini et al. 2021;Jelínek et al. 2022), but also to model various components of the GRB, such as forward shock (most of the works cited above employ MCMC), reverse shock (e.g., Laskar et al. 2016Laskar et al. , 2019) ) as well as prompt emission (Li 2022;Lazzati et al. 2023;Li et al. 2023).based on the optimisation of the search for suitable parameters sampling probability distribution functions (PDFs).
Parameters can be constrained by modelling the afterglow include the jet energy Ej, the jet opening angle θj, and the density n0 of the circumstellar medium, in addition to microphysical parameters of the particle acceleration process, including the fraction of thermal energy in the population of accelerated electrons ϵe and in the magnetic field ϵB, and the power law index of the population of non-thermal electrons p (see, e.g., Sari et al. 1998;Granot & Sari 2002).In addition, the emission depends on the jet structure and the observer angle, specially for observers located off-axis (see, e.g., Salafia & Ghirlanda 2022, for a review).
Given the large number of parameters involved in the model, one of the main challenges present when fitting GRB data using models is the possible presence of degeneracy in the parameter estimation.The issue of degeneracy has been recognised as a potential challenge in several articles (Kim et al. 2017;Granot et al. 2018;Resmi et al. 2018;Troja et al. 2018;Wu & MacFadyen 2018;Fraija et al. 2019;Gill et al. 2019;Hajela et al. 2019;Laskar et al. 2019;Beniamini et al. 2020;Cunningham et al. 2020;Angulo-Valdez et al. 2023;Laskar et al. 2023;McDowell & MacFadyen 2023), but its effects on the predictions of MCMC methods have not been discussed in detail.Parameter degeneracy occurs when different sets of parameters reproduce the same set of system's behaviour.This can happen not only in the obvious case in which there are more parameters than observations, but also in the more subtle (and often difficult to identify) case in which several observations provide the same information, all together.
In this article, we examine both degenerate and nondegenerate models to assess the reliability of results obtained through the use of the MCMC method.Additionally, we apply it to the analysis of the GRB 170817A afterglow emission, as a special case where the degeneracy is highly evident.We will show that, if not considered carefully, parameter degeneracy can go undetected when applying MCMC techniques, and lead to potential misleading interpretations of the results.We stress, anyway, that the issue of parameter degeneracy is not limited to the MCMC algorithm or GRB 170817A itself, and it is common to most optimisation methods.Therefore, the same procedure and discussion presented here can be extended to any area where the number of parameters in the model is larger than the number of observables, which subsequently leads to the emergence of degeneracies.
Our paper is organised as follows.In Section 2, we present a model of the afterglow phase of GRBs.In Section 3, we analyse in detail the results obtained by the application of MCMC to GRB 170817A.Finally, in Section 4, we discuss the advantages and drawbacks of the MCMC method, the implications for the analysis of GRB 170817A, and summarise our conclusions.

Analytical Model
In this section, we present the result of fitting degenerate models using MCMC.We consider as an example a simplified model of the GRB emission during the afterglow phase.In the next section, we show that the same problem is present when interpreting observations of the GRB170817A, among others.
GRB afterglows are typically modelled by considering synchrotron radiation from a population of electrons accelerated by shocks (e.g.Miceli & Nava 2022).The dynamics of a decelerating spherical blast wave is fully determined by considering its isotropic energy (E0) and the density (n) of the environment.The synchrotron emission is typically described as a function of the microphysical parameters ϵe, ϵB and p.These parameters give the fraction of thermal energy, magnetic energy and the power-law index of the electron energy distribution for the shock.
In general, the emission also depends on the characteristics of the emitting region, i.e. the jet opening angle (θj) in the case of a top-hat jet, or the full jet structure (specified by the energy E(θ) or shock Lorentz factor Γ sh (θ)) in the case of a structured jet.Here, θ represents the polar angle, which is measured from the direction of the main axis of the jet.The function E(θ) is usually defined (in analytical models) by considering a specific angular dependence, that is, top-hat (where E(θ) = E0 when θ < θj and zero otherwise, being θj the jet opening angle), Gaussian, or power-law jet structures (see e.g.Kumar & Granot 2003;Beniamini & Nakar 2019;Beniamini et al. 2020), among others, or calculated directly from the results of numerical simulations (see e.g.Lazzati et al. 2018;Xie et al. 2018;Salafia et al. 2020;Gottlieb et al. 2021;Urrutia et al. 2023).Setting the jet structure substantially reduces the number of independent parameters.Furthermore, for jets observed off-axis (as the GRB 170817A), or when modelling observations made after the jet break time, it becomes necessary to also consider the position of the observer, specified by the observer angle (θ obs ).
The normalisation, slope, and general behaviour of the spectrum at a particular time depend on the characteristic frequencies.Inspired by observations of GRB 170817A, we consider a spectrum that lies entirely within the same spectral range at the specific time considered.In particular, we consider the case νm ≪ ν ≪ νc, where the frequencies νm and νc correspond to the minimum and maximum energy of the electron population (e.g., Sari et al. 1998;Granot & Sari 2002).The time evolution of the flux at a given frequency depends on the details of the jet dynamics and the observer angle.For the dynamics described by the Blandford & Mc-Kee (1976) self-similar evolution for an on-axis observer in an homogeneous environment, in particular, we have (Granot & Sari 2002): Following Granot & Sari (2002), we indicate this spectral range as the PLS (power-law segment) "G".
To see the impact of employing an MCMC fitting procedure on a degenerate problem, we further simplify equation 1, by writing it as Then, we set and we fit synthetic data generated using equation 2. We set the constant 1.51 × 10 6 to obtain typical optical fluxes, corresponding to AB magnitudes between 16 and 22 (see e.g.Becerra et al. 2023b), and to explicitly show that k and ϵe are degenerate.
Figure 1 shows the synthetic spectrum generated by sampling data across radio frequencies (ν = 1.43 − 35 GHz in our example), optical (ν ∼ 1 eV), and X-ray frequencies (ν ∼ 1 KeV), and adding Gaussian noise at a 10% level.It is evident that, given the observed (synthetic) data points and the model described by the equation 2, the parameters k and ϵe remain strongly degenerate.In fact, while the equation 2 depends on three parameters (in addition to frequency), only two independent parameters can be determined from the synthetic spectrum, i.e. the slope of the spectrum (which implies that p will be well constrained by the fitting procedure) and the flux normalisation f .

Implementation of the MCMC method
We perform the MCMC fitting using the Python package emcee (Foreman-Mackey et al. 2013, 2019) and the guidelines from Appendix A.
To check that the MCMC method provides the correct result for a non-degenerate model, we first fit our synthetic data points by using the non-degenerated model described by: where the parameters to be fitted correspond to a constant k and p.
We impose the conditions 2 ⩽ p ⩽ 3 (Achterberg et al. 2001;Shen et al. 2006;Levan 2018) and k > 0 as priors for the log-prior function.We present the best-fit results (see Figure 1), the corner plot obtained for the MCMC procedure (see the top panel of Figure 2), and the autocorrelation plot k = ( 1.65 +0.17 0.16 ) × 10 6 (see the bottom panel of Figure 2).The autocorrelation time provides an estimate of the maximum number of initial iterations that should be discarded before the system reaches equilibrium (Sokal 1997;Foreman-Mackey et al. 2013).
As expected, the results exhibit robust convergence, with the best fit values closely approximating the original parameters used to generate the synthetic data points, consistently falling within a confidence interval 1σ.The value of k is determined with an error of 10%, consistent with the level of noise imposed on the data, while the value of p has a much smaller uncertainty (∼ 0.3%), due to the wide frequency range considered.
As can be seen in the bottom panel of Figure 2, the autocorrelation values for both p and k, followed over a span of 14200 iterations, fall below the 100N threshold imposed (see Appendix A).Therefore, we can state that, at the end of the iterative process, the Markov Chains related to the MCMC routine are independent of the initial position selected in the optimisation process, and the best-fit values obtained do not depend on the initial conditions.
Analogously, Figure 3 presents the results of the MCMC fitting procedure obtained in the degenerate case (i.e., with the model specified by equation 2).As in Figure 2, we present here the corner plot (top panel) and the autocorrelation values estimated for each parameter and their average (bottom panel).For the priors, we impose the conditions 2 ⩽ p ⩽ 3 (e.g., Achterberg et al. 2001;Shen et al. 2006;Levan 2018), 10 −5 ⩽ ϵe ⩽ 0.3 (Beniamini & van der Horst 2017), and 3.5 × 10 9 ⩽ k ⩽ 6.5 × 10 92 .We use a logarithmic sampling, although similar results are obtained when using a linear sampling.
The MCMC procedure satisfies the autocorrelation criteria after 10600 steps (see the bottom panel of Figure 3).The corner plot shows that the slope p is well constrained (p = 2.165 ± 0.006), and actually identical to the one shown in Figure 2.Although k is poorly determined and the distribution of its value does not show a clear peak in the contour levels (see the top panel of Figure 3), ϵe is well constrained (ϵe = (1.1 ± 0.2) × 10 −3 ).This result is clearly incorrect and would lead, in a real physical problem, to a wrong estimation of the value of this parameter.Due to condition 3.5 × 10 9 ⩽ k ⩽ 6.5 × 10 9 and the degeneracy of the model, all values 0.00074 < ϵe < 0.00126 are equally probable, as long as the product kϵ p−1 e remains constant.Therefore, the distribution of the values of ϵe and k should be flat in the range considered.An example of the correct flat distribution that should be obtained in a degenerate problem is shown, for example, in the middle panel of Figure B1.
In summary, achieving a well-defined "peaked" distribution for a specific variable does not necessarily guarantee the precision of the estimation of this particular parameter, at least when dealing with degenerated systems.When modelling simple physical systems, choosing a range of prior values too narrow can result in erroneous parameter estimations.In the next section, we show that the same problem is present when considering more complex degenerate problems, regardless of the width of the selected priors.

FITTING THE GRB 170817A
GRB 170817A, a low-luminosity short GRB, was detected by the Fermi and INTEGRAL telescopes beginning 1.7 s after the gravitational wave (GW) signal observed by the LIGO and Virgo interferometers on 17 August 2017 (Abbott et al. 2017a,b).Both events are the result of the merger of two neutron stars in NGC 4993, a lenticular galaxy located at a distance of 40 Mpc.The afterglow emission from the GRB 170817A can be successfully interpreted as the emission of a structured jet observed off-axis (see, e.g.D 'Avanzo et al. 2018;Margutti et al. 2018;Troja et al. 2018).Makhathini et al. (2021) presented the full multiwavelength afterglow light-curve data of GRB 170817A between 0.5 and 940 days after the merger, compiling all available photome- try3 .Makhathini et al. (2021) fit the available data using a smoothly broken power-law model The six parameters of their model (Fp, s, α1, α2, tp and β) are determined precisely by using MCMC (see their figure 3).Consistent with this phenomenological model, the GRB 170817A spectrum falls, at all times, in the PLS G considered in the previous section (Margutti et al. 2017;Troja et al. 2019).Thus, the parameter β can be taken as constant as a function of time.Indeed, it is easy to see that just six parameters can be determined from the observations.These correspond to the slope of the spectrum β, the slope of the light curve as a function of time before and after the peak (α1 and α2 in equation 5), the flux normalisation Fp, the time tp corresponding to the peak in the flux, and the slope s of the light curve close to its maximum.As this model fits well the available data, other models which include extra parameters (e.g., including a more complex description of the rise or fall of the light curve) will lead to poorly constrained or degenerated parameters.Furthermore, flux normalisation depends on ϵe, ϵB, E0, and n.Thus, we expect that there will be a complete degeneracy between these parameters.This is illustrated in detail in the following.
Figure 4 shows our fit of the photometry observations of GRB 170817A with the afterglowpy library (Ryan et al. 2020).Figure 5 shows the corner plot of our fit of the GRB 170817A afterglow, and Figure 6 shows the evolution of the autocorrelation function for all parameters and their average.The values of the parameters inferred from the fit are shown in Table 1.The model depends on the isotropic energy of the jet E, the density of the circumstellar medium n (computed assuming a uniform medium), the microphysical parameters associated with the particle acceleration process ϵe, ϵB, p, the observer angle θ obs , and parameters which determine the jet structure.We consider a Gaussian jet model, which depends on θcore and θwing, respectively, the angle at which most of the jet energy is contained and the maximum angle of the jet.
As mentioned above, we fix p = 2.168 following Makhathini et al. (2021).As the value of p is very well determined by spectral observations ranging from radio to X-ray frequencies, fixing it has the only effect of making the convergence of the fitting process faster and does not affect the results.We also fix the fraction of accelerated electrons χ = 1 to limit the number of parameters.The priors were selected for a wide ranges of values, i.The model fits the data satisfactorily (see Figure 4).The corner plot illustrates whereas the different parameters are well constrained and how they are correlated with each other (by looking at the covariances between the different parameters).Figure 5 shows that some of the parameters are poorly constrained.For instance, the distribution of θwing is nearly flat in the range (0.4, 1).This is not surprising, as the energy at the edge of a Gaussian jet is very low.Thus, the emission from the jet at angles θ ≲ θwing is weak and does not substantially modify the light curve.Although noisy, the covariance plots show that θwing is not independent of other parameters.We can conclude that there is some level of degeneracy between θwing and the other parameters, and the estimation of θwing = 0.65 +0.24 −0.27 is not reliable.We now focus our attention on other parameters, which seem well determined by looking at the corner plot.Nakar & Piran (2021) have discussed in detail the degeneracy between θ obs and θcore, showing that different combinations of these parameters, together with the shock Lorentz factor (which depends on the isotropic energy and ambient density here) produce the same results.Figure 5 shows that there is a certain degree of correlation between the values of θj and θcore, although they are both determined with a small error and do not appear to be fully degenerate from the figure.
The projection of the posterior probability distributions of energy and density also presents a clear peak.Furthermore, the covariance between these two parameters clearly restricts the range of possible values to a small region of the parameter space (around E0 = 10 52.53 erg and n = 3.1 × 10 −3 cm −3 ).Thus, we conclude by looking at the fit and the corner plot that density and energy are also well constrained and not degenerated.
The results shown in the figure 4 and discussed in the last two paragraphs are both incorrect and can lead to a wrong interpretation of the physical conditions of this GRB.To verify the degeneracy between energy and density, in particular, we show in Figure 4 the afterglow light curve produced by three afterglowpy models with different values of E0, n0 and ϵB.Indeed, the light curves look nearly identical.This is consistent with the scaling expected in the particular spectral range in which the GRB 170817A observations fall at all times.The peak flux and peak time, in fact, scale as (van Eerten & MacFadyen 2012; Granot 2012): That is, any combination of E, n, ϵe and ϵB that maintains unchanged tp and Fp, will produce the same light curve.This can be seen in Figure 4, in which we have multiplied/divided  5.
the best value of energy and density by the same amount (a factor of three in this case) between the different models.
As energy and density are rescaled by the same amount, equation 7 shows that, to maintain the same peak in flux, the scaling in energy must be related to the scaling in the microphysical parameters by the equation To get the parameters shown in Figure 4, thus, we fix the value of ϵe, we increase/drop ϵB by a factor of 10, and compute using this equation the resulting scaling in energy4 .The scaling in density is taken equal to the one in energy, to maintain the position of the peak constant5 .Figure 6 shows the autocorrelation plot.As we can see, the autocorrelation plot does not show convergence either after 100 thousand steps.In this case, this plot implies that some of the parameters can be degenerated, and the results should be carefully considered.
In conclusion, GRB 170817A, one of the most recently studied events, is an example of how a model of a light curve of a GRB using MCMC could lead to misleading results.Despite the large number of observations, the determination of the density, energy, and microphysical parameters for this event, as well as other parameters, remains uncertain.

DISCUSSION AND SUMMARY
In the previous section, we have shown that inferring parameters from a degenerate model leads to misleading results, which can go undetected by the MCMC algorithm.Although a detailed analysis of the MCMC algorithm is outside the scope of this paper, as we do not pretend to present a full analysis of the advantages or drawbacks of this method, we have shown that: a) the MCMC method correctly gives the correct solution in non-degenerate and simple degenerate problems.b) When at least some of the priors are defined over a small range, degeneracy on other parameters can go unrecognised by the MCMC method (see section 2.2).In this case, the MCMC may converge to a specific solution, instead of recognising the degeneracy present in the model.c) The autocorrelation can be used as a criterion to determine whether a model is adequate or not (Roy 2020;Brooks & Roberts 1998), as long as the range of priors is not too narrow.However, as shown in section 2.2, the opposite is not true, and the convergence of the autocorrelation does not necessarily imply that the results converge to the correct solution.
Therefore, in degenerated systems, a detailed analysis of the model is needed to understand which (if any) variable is degenerated, and to understand how many parameters can be well constrained by the data itself.

How to break the degeneracy: the case of GRB 170817A
Given the number of parameters involved in models used to explain the phenomenology of GRBs, it is expected that in most cases there will be some degeneracy between the model parameters.Without more observational data than photometry, it is very difficult to restrict the parameters of the model that describe the dynamics and emission properties of the system.
Although fixing some of the parameters does not break the degeneracy, it helps to understand the relationship between the different parameters in the model.This is what has typically been done (implicitly) when comparing observations with numerical simulations.Relying on a set of underlying hypotheses on the energy, structure of the jet launched from the central engine, density of the dense environment (the star in the case of long GRBs or the debris of the neutron star merger for short GRBs), and so on, some of the parameters can be directly inferred from the result of the numerical simulation (e.g., the jet structure), and the observations can be then fitted, reducing the degree of degeneracy (e.g., Margutti et al. 2018;Wu & MacFadyen 2018;Xie et al. 2018, in the case of GRB 170817A).
Although the problem of parameter degeneracy is often present when modelling GRB data, its relevance in limiting the parameters of the model is often underestimated when fitting data by using MCMC.While a degenerate model should have as a solution a flat distribution of probability, we have shown that parameter degeneracy may not appear clearly in corner plots usually used to interpret the model fitting process.This is the reason why, when the problem of degeneracy has been discussed in detail, at least in the context of GRB 170817A, this has been done mainly by starting from a theoretical understanding of the model, instead than starting from the result of the MCMC process.
Nakar & Piran (2021) discussed the degeneracy between θcore and θ obs , and, in general, the jet structure.They have reported different estimations present in the literature regarding these two parameters and have explained the differences as due to the different choice of priors.As we have shown, the convergence in the case of the afterglow fit of GRB 170817A can be extremely slow, which implies that, in fact, the final result can still depend on the initial choice of the parameters.This can be potentially verified by plotting the autocorrelation plot together with the corner plot.
GRBs seen on-axis typically provide limited information on the jet structure or on the observer angle, as most of the emission comes from the core of the jet, and it is weakly dependent on the emission coming from larger polar angles (but see, e.g., Cunningham et al. 2020;Gill & Granot 2023;O'Connor et al. 2023 for cases in which the jet structure is determined from on-axis observations).On the other hand, degeneracy in the model is particularly problematic in GRBs seen off-axis as the model will also depend on the jet structure.In the case of GRB 170817A, the full spectrum falls in the same spectral range (PLS G, corresponding to νm ≪ ν ≪ νc), which implies that some of the parameters related with νm and νc (ϵe, ϵ b , E0, and n0, see Granot & Sari 2002) have a low value.Getting a spectrum that covers different power-law segments of the spectrum (e.g., Harrison et al. 2001;Chandra et al. 2008;Perley et al. 2014;Tanvir et al. 2018;Salafia et al. 2022;Wang et al. 2022) can substantially, if not completely, reduce the degeneracy 6 .
GRB 170817A is one of the most studied GRBs.Makhathini et al. (2021) presented the full data set, modelling the data with a phenomenological model.Their model includes six parameters, equal to the number of observables which can be inferred from the data set.Actually, most of the papers on GRB 170817A fitted more than six parameters.As shown in Section 3, there is a full degeneracy between ϵe, ϵB, E and n.Several authors (Granot et al. 2018;Wu & MacFadyen 2018;Gill et al. 2019;McDowell & MacFadyen 2023) have discussed in detail the high level of degeneracy present in the model, considering the same scaling between the physical quantities discussed in this paper.Granot et al. (2018) concluded that the number of parameters outnumber the number of constraints in the data, and a significant level of degeneracy also remains when rich multiwavelength data are available, as in the case of GRB 170817A.Gill et al. (2019) also discussed in detail the degeneracy related to the scaling of the light curve.They determined a lower limit on the jet energy of 5.3 × 10 48 ergs, and on the ambient density of 5.3 × 10 −6 cm −3 , also considering observations of Very Long Baseline Interferometry (VLBI).By considering the diffuse X-ray emission from the hot plasma in the host galaxy, (Hajela et al. 2019) estimated an upper limit on the ambient density of 9.6 × 10 −3 cm −3 .Using this estimation for the ambient density, McDowell & MacFadyen (2023) estimated a jet energy of 7.5 × 10 48 erg.Using the value of the density estimated by (Hajela et al. 2019) as an upper limit, instead, we get the following ranges for the physical parameters: E = 8 × 10 50 − 1 × 10 53 erg, n = 7 × 10 −5 − 9 × 10 −3 , ϵe = 10 −3 − 0.3, ϵB = 10 −10 − 0.3.
The degeneracy can also be reduced when observations of polarisation are available, or the jet is resolved spatially, i.e. by VLBI observations.The upper limit of 12% on the linear polarisation of GRB 170817A available at 244 days in radio frequencies was used to constrain the magnetic field anisotropy factor B ∥ /B ⊥ = 0.5 − 0.9 both from analytical methods and numerical simulations (Gill & Granot 2020;Medina Covarrubias et al. 2023), consistent with previous estimates (Granot & Königl 2003;Stringer & Lazzati 2020).In principle, polarisation can also place constraints on the other parameters of the afterglow model.The same is true if VLBI observations are available (see, e.g., Taylor et al. 2004;Mooley et al. 2018a;Ghirlanda et al. 2019;Salafia et al. 2022).
Another approach to reduce the degeneracy is to use a smaller subset of parameters that can be constrained without a global fit of the GRB light curves (see, e.g.Santana et al. 2014;Beniamini & van der Horst 2017).For instance, as demonstrated by Beniamini & van der Horst (2017), it is possible to constrain ϵe by using the peaks of radio light curves.
Finally, we note that, in the era of multi-messenger astrophysics, simultaneous measurements of GRBs with gravitational waves and/or neutrinos allow search intervals to be eventually reduced (e.g.see Finstad et al. 2018;Lyman et al. 2018;Dichiara et al. 2021).

Final thoughts
In this paper, we have presented several examples to understand the effect of fitting the afterglow of GRBs using the MCMC method.The ease of implementation, as well as the availability of different free libraries, has made this method very popular in the analysis of astrophysical data.Despite the widespread use of this technique, the degeneracy eventually present in the model has not been thoroughly considered in the GRB literature.
In this work, we have emphasised the importance of considering the possible parameter degeneracy and implementing appropriate priors when interpreting the results obtained by employing the MCMC algorithm.Although the computational time required for correct modelling observations can be considerable (South et al. 2022), we highlighted the importance of conducting a convergence analysis, which is necessary to demonstrate that the derived interpretations are reliable (Roy 2020;Brooks & Roberts 1998).In cases of parameter degeneracy, the lack of convergence can be clearly observed in the autocorrelation plot (see Figure 6).This occurs because the MCMC can get trapped in a region of high probability density that may not necessarily represent the correct physical model.It is also important to note that MCMC can converge if the initial range of priors has been chosen too narrow.In such cases, MCMC may find a solution that fits the data well, but in reality is constrained by the initial priors and may not accurately represent the correct physical model.
Regarding GRB 170817A, we have shown that it is not possible to determine the energy of the explosion, the density of the environment and the microphysical parameters ϵe and ϵB related to the particle acceleration process from the afterglow light curve alone.Furthermore, for the specific case of GRB 170817A, the complete radio to X-ray afterglow spectrum lies on the same spectral segment making the degeneracy problem much more severe.Nevertheless, the degeneracy can be identified in most GRBs.a few examples include GRB 230812B (Hussenot-Desenonges et al. 2023), GRB 221009A (Laskar et al. 2023), GRB 210205A (Gupta et al. 2022), GRB 210104A (Zhang et al. 2022), or GRB 160625B (Cunningham et al. 2020).
Dealing with parameter degeneracy is very difficult not only for the MCMC method but more generally for most optimisation methods.This highlights the importance of fully understanding the possible presence of parameter degeneracy in a model during the fitting process.There are several practical steps that can be taken to avoid misinterpreting the results of MCMC, including using multiple Markov chains using different priors7 , implement other diagnosis tools such as the Gelman and Rubin statistics (Gelman & Rubin 1992), or by choosing priors that better reflect previous knowledge or assumptions about the parameters, rather than using flat priors.Moreover, it is recommended to use alternative fitting algorithms such as local and global minimizers, nested sampling or machine learning methods (Kochoska et al. 2020).
Finally, the degeneracy problem extends beyond the analysis of GRBs, and the discussion presented in this work is applicable to any other field in astrophysics that involves model fitting with a limited number of observables.In such cases, it is recommended to exercise particular caution with employing fitting techniques like MCMC, specially in situations where parameter degeneracy may remain unnoticed.Failing to detect degeneracy might result in misleading results, analysis and interpretations.
assumptions or information about the parameters, and it imposes constraints on the parameters of a model (MacKay 2002).The log-likelihood function calculates the acceptance probability for each step within the parameter space for every chain.(i.e., it measures the accuracy of fitting a statistical model), and the log-probability function maps the parameter space and is used to calculate the PDF in bayesian inference (MacKay 2002;Foreman-Mackey et al. 2013).
In this work, to avoid making assumptions on the a priori distribution, we use flat priors on all parameters in every MCMC routine.Moreover, we use a like-chi square formula as log-likelihood function: where ti and νi refer to the time and frequency data i, respectively, Fe,i is the uncertainty associated with each value of the observed flux Fi, F (θ, ti, νi) corresponds to the model function to be fitted to the data, and θ are the set of parameters associated (see section 2.1).
In MCMC approaches, one of the most important parameters is the number of workers in the process.This parameter represents the number of parallel chains running simultaneously in the MCMC routine, which are used to explore the parameter space on the sampler (Foreman-Mackey et al. 2013).In general, it is recommended to use hundreds of workers, although it is possible to go beyond (i.e., thousands of workers).In this paper we use 1000 workers to obtain a sufficient number of independent samples per autocorrelation time (whose value provides an estimate of the maximum number of initial iterations that should be discarded before the system reaches equilibrium (Sokal 1997;Foreman-Mackey et al. 2013).
A possible misuse of emcee is to set a process with a fixed (and random) number of samples.As Foreman-Mackey et al. ( 2013) recommended, we only need to run the sampler a few autocorrelation times.Beyond that limit, we obtain an independent set of samples.With that in mind, we monitor the progress of the chain using an autocorrelation time criterion τ that is updated every 100 steps (Foreman-Mackey et al. 2013).Therefore, we stop the sampler if the Markov chain has a length L with L ⩾ 100τ and if τ has not changed by more than 1% with respect to the previous estimation of τ (100 iterations before).In this way, the convergence of the Markov chain is ensured, removing the need to use a fixed number of steps.However, we set a maximum number of iterations of 10 5 to handle cases where there is no clear convergence (i.e., degenerate models with a very low convergence rate).
By employing an optimisation step before sampling can substantially mitigate the computational cost associated with posterior probability density functions, leading to more enhanced and well-suited fittings.Therefore, we use a maximum likelihood estimation before executing the MCMC sampling process.
Moreover, it is important to note that MCMC is a sampler, not an optimiser (Hogg & Foreman-Mackey 2018).Therefore, if an optimisation of the posterior probability density functions (PDFs) or likelihood is required, another process will be needed before implementing MCMC.

Figure 1 .
Figure 1.Best fit (red line) obtained using MCMC with the model function specified by equation 1. Synthetic data points (black points) are generated using the same model, adding a Gaussian with Gaussian noise at a 10% level.

Figure 2 .
Figure 2. Top panel: Corner plot obtained in MCMC on the nondegenerated model function (equation 4) using emcee (Foreman-Mackey et al. 2013).Labels k, p correspond to the fitted coefficients in the model.The values used to generate the synthetic data are shown with red lines.Bottom panel: Autocorrelation plot obtained for the MCMC routine of the model in equation 2. The black dashed line represents the threshold N = 100n, with n the number of iterations of MCMC.The red line represents the autocorrelation time τ estimated by the emcee, and the dotted lines represent the autocorrelation time estimated on the parameters k (blue line) and p (green line) as subindices of τ .

2Figure 3 .
Figure 3. Top panel: Corner plot obtained in MCMC on the degenerated model function (equation 2) using emcee (Foreman-Mackey et al. 2013).Labels k, ϵe, p correspond to the fitted coefficients in the model of Eq. 2. The original values provided for the synthetic data are shown with red lines.Bottom panel: Autocorrelation plot obtained for the MCMC routine of the model in equation 2. The black dashed line represents the threshold N = 100n, with n the number of iterations of MCMC.The red line represents the average autocorrelation time τ estimated by emcee.The dotted lines represent the autocorrelation time for the parameters α (blue line), ϵe (green line) and p (magenta line).

Figure 4 .
Figure 4. Fit of the GRB 170817A afterglow.The data points represent the observations, rescaled at ν obs = 3 × 10 9 Hz.The models are obtained for different combinations of the isotropic energy of the jet E, the density n of the environment, and the microphysical parameter ϵ E , and illustrate the intrinsic degeneracy in the fitting process.Data taken from Makhathini et al. (2021).

Figure 5 .
Figure 5. Corner plot of our fit of GRB 170817A afterglow.The values of the best fit obtained for each parameter are coloured red and listed in Table 1.Data were taken from Makhathini et al. (2021).

Figure 6 .
Figure6.Autocorrelation plot of our fit of GRB 170817A afterglow showed in Figure5.This figure illustrates that there is no convergence for any of the parameters.

Figure B1 .
Figure B1.Top panel: Synthetic data points (black lines) and best fit (red line) using MCMC with the model function of equation B2 of the data with Gaussian noise (black points).Central panel: Corner plot obtained in MCMC on the degenerated model function (equation B2) using emcee (Foreman-Mackey et al. 2013).The labels a, b, c correspond to the fitted coefficients in the model.The true values used to generate the synthetic data are shown with red lines.Bottom panel: Autocorrelation plot obtained for the MCMC routine of the model in equation 2. The black dashed line represents the threshold N = 100n, with n the number of MCMC iterations.The red line represents the average autocorrelation time τ estimated by emcee.The magenta dotted line represent τ for parameters b and c, and the blue dash-dotted line represent τ for parameter a.

Table 1 .
Parameters of the best fit model of GRB 170817A shown in Figure