From photons to neutrinos : Advancing multi-messenger modeling with Gammapy

Introduction

The study of extreme particle accelerators requires collecting and analyzing data from a large number of instruments, covering a large part of the photon spectrum, or combining photons with additional messengers such as neutrinos or gravitational waves. Current approaches for multi-wavelength and multi-messenger analyses typically involve the computation of flux points estimated in each waveband, followed by fitting a physical model with $\chi^2$ statistics. This method has, however, several limitations. First, flux points are derived with the assumption that the flux in each energy band follows a functional form (such as a Power Law or a Log Parabola), which is different from the model that will later be fitted to these flux points, introducing a bias. Additionally, upper limits (ULs) are often only used as strict thresholds, to ensure that the model does not overshoot them.

The software Gammapy [1] is the core library for the high-level analysis of CTAO data. Its extension, gammapy-mwl¹, enables the analysis of data ranging from optical to gamma-rays. We built a framework based on this extension with the version 2.0 of Gammapy [2] that allows for statistically sound multi-wavelength fitting of arbitrary models through forward-folding to instrumental counts. An approach using functional forms fitting on multi-wavelength data has been presented in [3]. Because this framework operates directly on events/counts, it inherently accounts for ULs and does not introduce any bias from flux-point estimation. The performance of this tool is presented in the following sections.

An application: simulating blazar emission

We created two types of datasets to probe the performance of a fit to counts compared to a standard fit to flux points. We simulated the count datasets with real instrument response functions (IRFs) of UVOT, XRT and LAT and the following simulated IRF of CTAO Prod5-North-40deg-AverageAz-4LSTs09MSTs.18000s-v0.1.fits[4]. The assumed spectral model was a Synchrotron Self-Compton (SSC) model of a blazar, simulated using the Gammapy plugin of the JetSeT code [5, 6, 7]. The model includes optical reddening, absorption by neutral hydrogen in soft-X-rays and absorption by pair-production on the extragalactic background light. As a first test, we simulated the spectral energy distribution of a typical high- to extreme-high-frequency-peaked BL Lac object, using model parameters close to the ones used in the modeling of the well-known blazar Mrk 501 [8]. The simulated count datasets were created from UVOT and XRT data from a single Swift pointing of 450 s, Fermi/LAT data integrated over 10 days, and a 1-hour simulation of CTAO observations (using the North array). We generated flux points for each instrument separately, and merged them into a flux-point Gammapy dataset, as it is the current standard in multi-wavelength studies. Within each individual band, we adopted as spectral model a power-law function, or a log parabola if it was favored by more than 3 $\sigma$ over the power law. An example of the simulated count dataset, and the flux-point dataset is shown in the left and right plots of Figure 1, respectively.

Performance of the two algorithms

We simulated 100 different pairs of count and flux-point datasets to study the average properties of the new framework. We then fit each dataset with the same starting parameters, using the SSC JetSeT plugin in Gammapy. As we expect convergence to be more difficult to achieve with a higher number of free parameters, we freed the SSC parameters in the following order: $B$, $N$, $\gamma_{\mathrm{break}}$, $p_1$, and $p$. The simulated values are $B = 0.1$ G, $N = 5$ cm$^{-3}$, $\gamma_{\mathrm{break}} = 5 \times10^5$, $p_1 = 3.5$, and $p = 2.5$. Furthermore, the count datasets possess additional degrees of freedom compared to the flux-point datasets because of free normalization of other sources in the Fermi field of view.

The minimization within Gammapy was done using iminuit [9]. Two approaches were tested for initializing the parameters:

• Good initial guess: Free parameters among those mentioned before are initialized with $B = 0.05$ G, $N = 10$ cm$^{-3}$, $\gamma_{\mathrm{break}} = 2 \times10^5$, $p_1 = 3.2$, and $p = 2.2$, while frozen parameters are set to their simulated values. This approach mimics a pre-fit of a model by eye, which we would typically attempt in order to estimate reasonable starting parameters before running a fit. In this case, all 100 fits for both flux-point and count datasets share the same starting parameters.

• Random initial guess: Starting parameters are drawn randomly following a uniform distribution. Here, each pair of flux-point and count datasets shares the same random starting parameters, but these differ across pairs.

Additionally, we derived posterior probability distributions assuming uniform priors and the Bayesian evidence with the nested sampling Monte Carlo algorithm MLFriends [10, 11] using the UltraNest² package [12]. We then compared the contours extracted from the covariance matrix, estimated with iminuit strategy 1, with the nested sampling posterior distributions. An example of a convergent fit of a single realization is shown in the left and right plots of Figure 2, for the count dataset and the flux-point dataset, respectively.

The fitted physical parameters yield very similar values whether we use the flux-point or count dataset. While the nested sampling posteriors align well with the covariance contours for the count dataset, the covariance contours for the flux-point dataset are significantly smaller than the nested sampling posteriors. This discrepancy suggests that errors are greatly underestimated when fitting to flux points. We can further support the former statement by looking at the number of realizations for which the fitted parameters are within a $5\sigma$ interval around the simulated parameters. We compute this interval through the Mahalanobis distance (1).

In this equation, $u$ and $v$ represent the vectors of fitted and simulated parameters, and $V$ is the covariance matrix.

Fitting to counts significantly increases convergence compared to fitting to flux points. Any conclusion on the quality of the fits is however non-trivial. To address the latter point, we evaluate how many of the 100 realizations produced convergent fits with reliable uncertainties and a Mahalanobis distance of less than 5 between the fitted and simulated parameters, which we call "accurate fits". The results of these two approaches are shown in Figure 3. The orange and blue bars respectively represent the count and flux-point datasets. Darker bars correspond to fits initialized with a good guess, while lighter bars represent fits initialized with random parameters. The reliability of

the fit to flux points is low even with a low number of free physical parameters. On the other hand, the fit to counts is more robust even with a high number of free parameters. In addition, starting with a good guess of the initial parameters matters, especially when several physical parameters are free.

In terms of computing time, the fit to counts is not significantly longer than the fit to flux points. Indeed, the standard fit to flux points with 3 free physical parameters takes $\sim5$ min, while the standard fit on counts takes $\sim10$ min. The nested sampling, which improves the uncertainty on the parameters for the fit to flux points, is much more computationally intensive, and requires over 3 h for both the flux point and count datasets.

Conclusion and outlook

We developed a novel framework within the open-source Gammapy software, dedicated to unbiased multi-wavelength and multi-messenger fitting. This framework enables fast and robust fitting of any physical model directly on instrumental counts. First applications show that it not only improves convergence but also provides statistically accurate errors on the physical parameters. Nested sampling, whether used on counts or flux points, remains the most accurate method to estimate parameter values and uncertainties, at the expense of a much longer computing time. However, multi-wavelength fitting with iminuit to counts provides a much faster convergence and solid estimate of uncertainties. The uncertainties on the fitted parameters estimated from the iminuit fit to flux points cannot be trusted. While our application focuses on a leptonic model, this tool can be used on (lepto-)hadronic models as well. Neutrino analysis is already compatible with Gammapy [3, 4] and will soon be officially implemented allowing for statistically sound multi-messenger fits. Since hadronic fits are computationally intensive, we are exploring the use of precomputed grids and machine learning solutions to reduce the computing time and make the fitting algorithm attractive for a large variety of astrophysical applications.

Acknowledgments

This work benefited from the support of the project COCOA-NuGETs ANR-23-CE31-0026 of the French National Research Agency (ANR). This work made use of Gammapy [1], a community-developed Python package. The Gammapy team acknowledges all Gammapy past and current contributors, as well as all contributors of the main Gammapy dependency libraries: NumPy, SciPy, Astropy, Astropy Regions, iminuit, Matplotlib.