# Detection of the Geminga pulsar at energies down to 20 GeV with the LST-1 of CTAO

K Abe, S Abe, A Abhishek, F Acero, A Aguasca-Cabot, I Agudo, C Alispach, D Ambrosino, F Ambrosino, L.A Antonelli, et al.

► To cite this version:

K Abe, S Abe, A Abhishek, F Acero, A Aguasca-Cabot, et al. Detection of the Geminga pulsar at energies down to 20 GeV with the LST-1 of CTAO. Astron. Astrophys., 2025, 698, pp.A283. 10.1051/0004-6361/202554350. hal-05108258

HAL Id: hal-05108258

https://hal.science/hal-05108258v1

Submitted on 16 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Detection of the Geminga pulsar at energies down to 20 GeV with
the LST-1 of CTAO

K. Abe¹, S. Abe², A. Abhishek³, F. Acero⁴,⁵, A. Aguasca-Cabot⁶, I. Agudo⁷, C. Alispach⁸, D. Ambrosino⁹, F. Ambrosino¹⁰, L. A. Antonelli¹⁰, C. Aramo⁹, A. Arbet-Engels¹¹, C. Arcaro¹², T. T. H. Arnesen¹³, K. Asano², P. Aubert¹⁴, A. Baktash¹⁵, M. Balbo⁸, A. Bamba¹⁶, A. Baquero Larriva¹⁷,¹⁸, U. Barres de Almeida¹⁹, J. A. Barrio¹⁷, L. Barrios Jiménez¹³, I. Batkovic¹², J. Baxter², J. Becerra González¹³, E. Bernardini¹², J. Bernete²⁰, A. Berti¹¹, I. Bezshyiko²¹, C. Bigongiari¹⁰, E. Bissaldi²², O. Blanch²³, G. Bonnoli²⁴, P. Bordas⁶, G. Borkowski²⁵, G. Brunelli²⁶,²⁷,*, A. Bulgarelli²⁶, M. Bunse²⁸, I. Burelli²⁹, L. Burmistrov²¹, M. Cardillo³⁰, S. Caroff¹⁴, A. Carosi¹⁰, R. Carraro¹⁰, M. S. Carrasco³¹, F. Cassol³¹, N. Castrejón³², D. Cerasole³³, G. Ceribella¹¹,*, A. Cerviño Cortínez¹⁷, Y. Chai¹¹, K. Cheng², A. Chiavassa³⁴,³⁵, M. Chikawa², G. Chon¹¹, L. Chytka³⁶, G. M. Cicciari³⁷,³⁸, A. Cifuentes²⁰, J. L. Contreras¹⁷, J. Cortina²⁰, H. Costantini³¹, M. Dalchenko²¹, P. Da Vela²⁶, F. Dazzi¹⁰, A. De Angelis¹², M. de Bony de Lavergne³⁹, R. Del Burgo⁹, C. Delgado²⁰, J. Delgado Mengual⁴⁰, M. Dellaiera¹⁴, D. della Volpe²¹, B. De Lotto²⁹, L. Del Peral³², R. de Menezes³⁴, G. De Palma²², C. Díaz²⁰, A. Di Piano²⁶, F. Di Pierro³⁴, R. Di Tria³³, L. Di Venere⁴¹, R. M. Dominik⁴², D. Dominis Prester⁴³, A. Donini¹⁰, D. Dore²³, D. Dorner⁴⁴, M. Doro¹², L. Eisenberger⁴⁴, D. Elsässer⁴², G. Emery³¹, J. Escudero⁷, V. Fallah Ramazani⁴⁵,⁴⁶, F. Ferrarotto⁴⁷, A. Fiasson¹⁴,⁴⁸, L. Foffano³⁰, S. Fröse⁴², Y. Fukazawa⁴⁹, S. Gallozzi¹⁰, R. Garcia López¹³, S. Garcia Soto²⁰, C. Gasbarra⁵⁰, D. Gasparrini⁵⁰, D. Geyer⁴², J. Giesbrecht Paiva¹⁹, N. Giglietto²², F. Giordano³³, N. Godinovic⁵¹, T. Gradetzke⁴², R. Grau²³, D. Green¹¹, J. Green¹¹, S. Gunji⁵², P. Günther⁴⁴, J. Hackfeld⁵³, D. Hadasch², A. Hahn¹¹, M. Hashizume⁴⁹, T. Hassan²⁰, K. Hayashi²,⁵⁴, L. Heckmann¹¹, M. Heller²¹, J. Herrera Llorente¹³, K. Hirotani², D. Hoffmann³¹, D. Horns¹⁵, J. Houles³¹, M. Hrabovsky³⁶, D. Hrupec⁵⁵, D. Hui², M. Iarlori⁵⁶, R. Imazawa⁴⁹, T. Inada², Y. Inome², S. Inoue²,⁵⁷, K. Ioka⁵⁸, M. Iori⁴⁷, T. Itokawa², A. Iuliano⁹, J. Jahanvi²⁹, I. Jimenez Martinez¹¹, J. Jimenez Quiles²³, I. Jorge Rodrigo²⁰, J. Jurysek⁵⁹, M. Kagaya²,⁵⁴, O. Kalashev⁶⁰, V. Karas⁶¹, H. Katagiri⁶², D. Kerszberg²³,⁶³, T. Kiyomot⁶⁴, Y. Kobayashi², K. Kohri⁶⁵, A. Kong², P. Kornecki⁷, H. Kubo², J. Kushida¹, B. Lacave²¹, M. Lainez¹⁷, G. Lamanna¹⁴, A. Lamastra¹⁰, L. Lemoigne¹⁴, M. Linhoff⁴², S. Lombardi¹⁰, F. Longo⁶⁶, R. López-Coto⁷, M. López-Moya¹⁷, A. López-Oramas¹³, S. Loporchio³³, A. Lorini³, J. Lozano Bahilo³², F. Lucarelli¹⁰, H. Luciani⁶⁶, P. L. Luque-Escamilla⁶⁷, P. Majumdar²,⁶⁸, M. Makariev⁶⁹, M. Mallamaci³⁷,³⁸, D. Mandat⁵⁹, M. Manganaro⁴³, D. K. Maniadakis¹⁰, G. Manicò³⁸, K. Mannheim⁴⁴, S. Marchesi²⁶,²⁷,⁷⁰, F. Marini¹², M. Mariotti¹², P. Marquez⁷¹, G. Marsella³⁷,³⁸, J. Martí⁶⁷, O. Martinez⁷², G. Martínez²⁰, M. Martínez²³, A. Mas-Aguilar¹⁷,*, M. Massa³, G. Maurin¹⁴, D. Mazin²,¹¹, J. Méndez-Gallego⁷, S. Menon¹⁰, E. Mestre Guillen⁷³, S. Micanovic⁴³, D. Miceli¹², T. Miener¹⁷, J. M. Miranda⁷², R. Mirzoyan¹¹, M. Mizote⁷⁴, T. Mizuno⁴⁹, M. Molero Gonzalez¹³, E. Molina¹³, T. Montaruli²¹, A. Moralejo²³, D. Morcuende⁷, A. Moreno Ramos⁷², A. Morselli⁵⁰, V. Moya¹⁷, H. Muraishi⁷⁵, S. Nagataki⁷⁶, T. Nakamori⁵², A. Neronov⁶⁰, D. Nieto Castaño¹⁷, M. Nievas Rosillo¹³, L. Nikolic³, K. Nishijima¹, K. Noda²,⁵⁷, D. Nosek⁷⁷, V. Novotny⁷⁷, S. Nozaki¹¹, M. Ohishi², Y. Ohtani², T. Oka⁷⁸, A. Okumura⁷⁹,⁸⁰, R. Orito⁸¹, L. Orsini³, J. Otero-Santos⁷, P. Ottanelli⁸², M. Palatiello¹⁰, G. Panebianco²⁶, D. Paneque¹¹, F. R. Pantaleo²², R. Paoletti³, J. M. Paredes⁶, M. Pech³⁶,⁵⁹, M. Pecimotika²³, M. Peresano¹¹, F. Pfeifle⁴⁴, E. Pietropaolo⁵⁶, M. Pihet⁶, G. Pirola¹¹, C. Plard¹⁴, F. Podobnik³, M. Polo²⁰, E. Prandini¹², M. Prouza⁵⁹, S. Rainò³³, R. Rando¹², W. Rhode⁴², M. Ribó⁶, V. Rizi⁵⁶, G. Rodriguez Fernandez⁵⁰, M. D. Rodríguez Frías³², P. Romano²⁴, A. Roy⁴⁹, A. Ruina¹², E. Ruiz-Velasco¹⁴, T. Saito², S. Sakurai², D. A. Sanchez¹⁴, H. Sano²,⁸³, T. Šarić⁵¹, Y. Sato⁸⁴, F. G. Saturni¹⁰, V. Savchenko⁶⁰, F. Schiavone³³, B. Schleicher⁴⁴, F. Schmuckermaier¹¹, J. L. Schubert⁴², F. Schussler³⁹, T. Schweizer¹¹, M. Seglar Arroyo²³, T. Siegert⁴⁴, G. Silvestri¹², A. Simongini¹⁰,⁸⁵, J. Sitarek²⁵, V. Sliusar⁸, A. Stamerra¹⁰, J. Strišković⁵⁵, M. Strzys², Y. Suda⁴⁹, A. Sunny¹⁰, H. Tajima⁷⁹, M. Takahashi⁷⁹, J. Takata², R. Takeishi², P. H. T. Tam², S. J. Tanaka⁸⁴, D. Tateishi⁶⁴, T. Tavernier⁵⁹, P. Temnikov⁶⁹, Y. Terada⁶⁴, K. Terauchi⁷⁸, T. Terzic⁴³, M. Teshima²,¹¹, M. Tluczykont¹⁵, F. Tokanai⁵², T. Tomura², D. F. Torres⁷³, F. Tramonti³, P. Travnicek⁵⁹, G. Tripodo³⁸, A. Tutone¹⁰, M. Vacula³⁶, J. van Scherpenberg¹¹, M. Vázquez Acosta¹³, S. Ventura³, S. Vercellone²⁴, G. Verna³, I. Viale¹², A. Vigliano²⁹, C. F. Vigorito³⁴,³⁵, E. Visentin³⁴,³⁵, V. Vitale⁵⁰, V. Voitsekhovskyi²¹, G. Voutsinas²¹, I. Vovk², T. Vuillaume¹⁴, R. Walter⁸, L. Wan², M. Will¹¹, J. Wójtowicz²⁵, T. Yamamoto⁷⁴, R. Yamazaki⁸⁴, Y. Yao¹, P. K. H. Yeung²,*, T. Yoshida⁶², T. Yoshikoshi², and W. Zhang⁷³

(Affiliations can be found after the references)

<footnote>* Corresponding authors: `lst-contact@cta-observatory.org`</footnote>
ABSTRACT

**Context.** Geminga is the third gamma-ray pulsar firmly detected by imaging atmospheric Cherenkov telescopes (IACTs) after the Crab and the Vela pulsars. Most of its emission is expected at tens of giga-electronvolts, and, out of the planned telescopes of the upcoming Cherenkov Telescope Array Observatory (CTAO), the Large-Sized Telescopes (LSTs) are the only ones with optimised sensitivity at these energies.

**Aims.** We aim to characterise the gamma-ray pulse shape and spectrum of Geminga as observed by the first LST (hereafter LST-1) of the Northern Array of CTAO. Furthermore, this study confirms the great performance and the improved energy threshold of the telescope, as low as 10 GeV for pulsar analysis, with respect to current-generation Cherenkov telescopes.

**Methods.** We analysed 60 hours of good-quality data taken by the LST-1 between December 2022 and March 2024 at zenith angles below 50°. Additionally, a new *Fermi*-LAT analysis of 16.6 years of data was carried out to extend the spectral analysis down to 100 MeV. Lastly, a detailed study of the systematic effects was performed.

**Results.** We report the detection of Geminga in the energy range between 20 and 65 GeV. Of the two peaks of the phaseogram, the second one, P2, is detected with a significance of $12.2\sigma$, while the first (P1) reaches a significance level of $2.6\sigma$. The best-fit model for the spectrum of P2 was found to be a power law with a spectral index of $\Gamma = (4.5 \pm 0.4)_{-0.6_{\sigma sys}}^{+0.2_{\sigma sys}}$, compatible with the previous results obtained by the MAGIC Collaboration. No evidence of curvature is found in the LST-1 energy range. The joint fit with *Fermi*-LAT data confirms a preference for a sub-exponential cut-off over a pure exponential, even though both models fail to reproduce the data above several tens of giga-electronvolts. The overall results presented in this paper prove that the LST-1 is an excellent telescope for the observation of pulsars, and improved sensitivity is expected to be achieved with the full CTAO Northern Array.

**Key words.** astroparticle physics – stars: neutron – pulsars: general – pulsars: individual: Geminga pulsar – gamma rays: stars

1. Introduction

Pulsars make up the majority of the Galactic gamma-ray sky, with almost 340 identified sources between rotation-powered and millisecond pulsars in the Third *Fermi* Large Area Telescope (LAT) catalogue of gamma-ray pulsars, 3PC (Smith et al. 2023). Their gamma-ray emission is usually explained with curvature radiation models, which are compatible with the double-peaked phase-folded light curves and the sub-exponential cut-offs at a few giga-electronvolts observed in the spectra of *Fermi*-LAT pulsars. Imaging atmospheric Cherenkov telescopes (IACTs) later detected emission above 25 GeV from the Crab pulsar (Aliu et al. 2008). Then, more surprisingly, the Crab was discovered to show pulsed emission above 100 GeV (VERITAS Collaboration 2011; Aleksic et al. 2012) and even up to 1.5 TeV (Ansoldi et al. 2016). More recently, H.E.S.S. Collaboration (2023) discovered pulses up to 20 TeV from Vela. The detection of tera-electronvolt pulsed emission challenges the capability of curvature radiation models in explaining the entire gamma-ray spectrum from megaelectronvolt to tera-electronvolt energies.

Despite the important achievement, current-generation IACTs have limited sensitivity for the detection of pulsars. The Cherenkov Telescope Array Observatory (CTAO) (Zanin et al. 2021) will be the next-generation facility for ground-based gamma-ray observations, consisting of two arrays, one in the Northern Hemisphere and one in the Southern, composed of telescopes of three different sizes. This configuration will allow it to cover a wide energy range from around 20 GeV to almost 300 TeV. In particular, Large-Sized Telescopes (LSTs) are designed to have the highest sensitivity in the energy range from 20 to 200 GeV (Abe et al. 2023), being ideal instruments for studying pulsars in the gap between *Fermi*-LAT and the tera-electronvolt domain, covered by other IACT arrays. The first LST of the CTAO-North array, dubbed LST-1, has already been built and is currently in the final commissioning phase.

Geminga (PSR J0633+1746) is the third pulsar detected above $5\sigma$ significance by IACTs<fnref n="1" /> (MAGIC Collaboration 2020). It is a well-known and studied source, considered the arch type of a middle-aged radio-quiet gamma-ray pulsar. It is also one of the closest pulsars to us, with an estimated parallax distance

of $d = 250_{-62}^{+120}$ pc (Faherty et al. 2007). The SAS-2 satellite discovered it as an isolated gamma-ray source in 1972 without any counterpart at other wavelengths (Fichtel et al. 1975). Only twenty years later it was possible to identify it as a pulsar thanks to combined results from X-ray observations by ROSAT (Halpern & Holt 1992) and gamma-ray data from EGRET (Bertsch et al. 1992). From these last results, the main rotation parameters were estimated to be a period of $P = 237$ ms and a period derivative of $\dot{P} \simeq 1.1 \cdot 10^{-14} \text{ s s}^{-1}$, from which they derived its characteristic age of $\tau_c \sim 300$ kyr and its spin-down power $\dot{E} \sim 3.5 \cdot 10^{34} \text{ erg s}^{-1}$. There is still no detection of Geminga at radio wavelengths, and the latest results report an upper limit on the average flux density of 0.4–4 mJy at a frequency of 111 MHz (Ershov 2021), where the range is due to different assumptions on the pulse duration. Despite being radio-quiet, Geminga is one of the few pulsars also detected at optical frequencies (Shearer & Golden 2002).

Geminga was detected within one year of operations by *Fermi*-LAT (Abdo et al. 2010), and, more surprisingly, emission above 15 GeV was reported by MAGIC Collaboration (2020). They detected the second emission peak of the pulsar at a significance level of $6.3\sigma$ and fitted its spectrum using a power law (PL) with spectral index $\Gamma = 5.62 \pm 0.54$ extending up to 75 GeV. This could be interpreted as the smooth transition from a regime dominated by curvature radiation to an inverse-Compton-dominated regime at higher energies.

In this paper, we report the detection of Geminga with the LST-1. Section 2 summarises the details of the data analysis; in Sect. 3 we present the main results obtained in this work, in Sect. 4 we outline the tests we performed to study systematic effects on the LST-1 analysis; finally, in Sect. 5 we discuss the results and present our conclusions.

2. Observations and data analysis

2.1. LST-1 observations and data analysis

LST-1 observations were performed in 20-minute-long runs in the so-called wobble mode (Fomin et al. 1994), with an offset angle of $0.4^\circ$ from the camera centre. The LST-1 observed Geminga as part of the commissioning program between December 2022 and March 2024 for almost 80 hours in total. For this analysis, we required observations in dark conditions

<footnote n="1">¹ Note that a fourth pulsar, PSR B1706-44, was reported with a $4.7\sigma$ significance by the H.E.S.S. Collaboration (Spir-Jacob et al. 2019).</footnote>
only, i.e. no dusk and down or moonlight, and performed a
quality selection of the data to reject the runs affected by tele-
scope problems or bad weather conditions. The method chosen
to assess the data quality is based on the analysis of the differen-
tial intensity spectra of the detected showers, which are mostly
induced by hadrons. The spectrum is characterised by a peak
followed by a power law decay. The method consists of fitting
the power-law region of the spectrum, then the selection is based
on a range of values derived empirically from a good sample of
observations. This approach is an evolution of the one used in
CTAO-LST Project (2023). The final sample was reduced to 60
hours of good-quality data taken at zenith distance (Zd) lower
than 50°, among which the majority (54 hours) were taken at
Zd < 25°, which helps to reduce the overall energy threshold in
the analysis.

The data were then reduced using version v0.10.7
of the cta-lstchain software library (López-Coto 2022;
López-Coto et al. 2024), following the procedure described in
CTAO-LST Project (2023). Three main characteristics of the pri-
mary particle that generates the shower were reconstructed: its
energy, incoming direction in the sky, and nature. The latter was
measured with the ‘gammaness’ parameter, an indicator of how
likely it is that the primary particle is a hadron (gammaness ~ 0)
or a gamma ray (gammaness ~ 1). The reconstruction is per-
formed using random forest methods trained on simulated Monte
Carlo (MC) data. The package lstmcpipe (Garcia et al. 2022;
Vuillaume et al. 2024), which uses lstchain, was employed to
process the simulated data. In particular, an MC production of
gamma rays and protons was simulated along a declination line
of 22.76° as seen from the LST-1 site, close to the declination
of Geminga in the sky. We also tuned the night sky background
in the simulations to match the level observed in the sky region
around Geminga. The instrument response functions (IRFs) of
the LST-1 were computed using an MC test sample simulated in
different sky positions along the declination line, and they were
subsequently interpolated with the telescope pointing to obtain
the correct IRFs for each observation run.

We adopted the source-dependent approach, for which the
source position is assumed a priori. In this way, we can achieve
a better performance in the reconstruction of low-energy show-
ers (CTAO-LST Project 2023). The main parameter introduced
when performing a source-dependent analysis is ‘alpha’, defined
as the angle between the main axis of the shower ellipse and the
line connecting the centroid of the image and the position of the
source.

To achieve a good agreement between the simulations and
the observed data, especially for the weaker signals, we applied
to both samples a cut on the image intensity of *I* > 50 photoelec-
trons (p.e.). This cut allowed us also to reject the lowest intensity
events, which are usually poorly reconstructed. To remove the
hadronic background, we also applied energy-dependent cuts on
alpha and gammaness based on a standard 70% selection effi-
ciency of the MC. Later, in Sect. 4, we study the systematic
effects on the analysis when varying these cuts.

The next step of the analysis chain was the computa-
tion of the rotational phase of the pulsar at which each
event was emitted. This was done using version v0.9.7 of
the PINT Python package (Luo et al. 2019, 2021) through the
PulsarTimingAnalysis software (Mas & Morcuende 2024).
The pulsar’s ephemeris were computed based on Fermi-LAT
data starting on 6 August 2008 and available online<fnref n="2" />. Finally, the
science products, i.e. the phaseogram and the spectral energy dis-

tribution (SED), were obtained using version v1.1 of Gammapy
(Donath et al. 2023; Aguasca-Cabot et al. 2023).

The phase regions adopted for this analysis are the ones
defined in MAGIC Collaboration (2020): P1 = [0.056; 0.161]
for the first peak, P2 = [0.550; 0.642] for the second, and OFF
= [0.700; 0.950] for the background phase region. We also con-
sidered the inter-peak (or Bridge) emission, defined as Bridge =
[0.161; 0.550].

## 2.2. Fermi-LAT observations and data analysis

Most of the gamma-ray emission from Geminga is concentrated below 10 GeV, where Fermi-LAT is more sensitive than the LST-1. To expand the spectrum over a broader energy range, we analysed 16.6 years of public Fermi-LAT data, starting on 4 August 2008 and lasting until 27 March 2025.

Fermi-LAT data were reduced and analysed using version v1.3.1 of the Fermipy package (Wood et al. 2017) with version v2.2.0 of the Fermi Science Tools and the P8R3_SOURCE_V3 IRFs. We chose a region of interest (ROI) with a side of 15° centred on the Geminga pulsar position (RA = 06h 33m 54.15s, Dec = +17° 46′ 12.91″) and binned in 0.1° bins, and the energy range 100 MeV–300 GeV. We filtered for good timing and nominal state of the detector. For energies below 10<sup>3.5</sup> ~ 3.2 GeV, we selected the events classified as SOURCE (event class 128) and type 32 (PSF3, the best quartile of the PSF distribution), with a maximum zenith distance of 90°. For energies above ~3.2 GeV, instead, we selected the SOURCE class events, but we included all event types 4, 8, 16, 32 (PSF0 + PSF1 + PSF2 + PSF3, corresponding to no filtering for the PSF quality), with a maximum zenith distance of 105°. This event selection allows us to limit systematics related to the PSF modelling and affecting a bright source such as Geminga at the lowest energies, while keeping sound statistics at the higher end of the spectrum.

The last step was to obtain the spectrum and the flux points
for the two peaks, for which we selected the same phase regions
as for the LST-1 analysis. We performed a binned likelihood
analysis in the range 100 MeV–300 GeV with eight logarithmi-
cally spaced energy bins per decade. Then, we used a spectral-
spatial model considering the sources reported in the 4FGL
catalogue (Abdollahi et al. 2020) and within a region of 20°.
We included in the fit the Galactic (gll_iem_v07.fits) and
isotropic background (iso_P8R3_SOURCE_V3_v1.txt) mod-
els, specifying the latter according to each different PSF class
separately. We first ran an optimisation to estimate the predicted
count rates of the sources in the sky region and froze all the
parameters to their catalogue values for sources with predicted
count values of less than one. We then modelled the Geminga
spectrum as a power law with sub-exponential cut-off (PLSEC),
using the PLSuperExpCutoff2 model, first with free normalisa-
tion, spectral index and exponential factor, but leaving the expo-
nential index at the 4FGL catalogue value. We also let free the
normalisation of the background sources within 10° and with at
least a 3σ significance, as well as the shape parameters of the
background sources within 7° with at least a 5σ significance,
and the spectral parameters of the second brightest source in the
ROI, 4FGL J0617.2+2234e. Starting from the maximum likeli-
hood values of the previous fit, we performed a further step by
freezing background sources, but freeing all Geminga param-
eters, including the exponential index. Finally, we obtained
the Geminga SED points by fitting a power law with spectral
index Γ = 2 in each energy bin, with the background sources
locked.

<footnote n="2">² https://www.mpp.mpg.de/~ceribell/geminga/index.php</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Phaseogram of the LST-1 observations of the Geminga shown over two rotational periods, with no cut in energy. The different phase regions (P1, P2, Bridge and background, or OFF) are highlighted in the plot. The average level of the background counts is reported as the horizontal dashed line. We also show the Li&Ma significance of both peaks and the inter-peak region and the total observation time.</figcaption></figure>

# 3. Results

## 3.1. Phaseogram

Figure 1 depicts the phaseogram of Geminga as observed by the LST-1 at $Zd < 50^{\circ}$ without applying any energy cut to the sample. Out of the two well-known peaks, only P2 is detected with a Li&Ma significance (Li & Ma 1983) of $12.2\sigma$. P1 reaches a significance level of $2.6\sigma$. We also studied the signal in the Bridge region and observed an excess corresponding to a Li&Ma significance of $1.8\sigma$.

We expect that the significance evolves as the square root of the observational time, i.e. $\sigma = A\sqrt{t(h)}$, so we used this relation to derive the coefficient A and estimate how many hours would be needed to reach a $5\sigma$ significance of P1. The test yielded $\sigma_{P1} = (0.34 \pm 0.13) t(h)^{1/2}$, where the error on A has been estimated considering that the Li&Ma significance, in the case of weak sources, has an uncertainty of one by definition (Li & Ma 1983). Considering this trend, the estimated observational time for a $5\sigma$ detection of P1 would be more than 200 hours, which could be achieved with few observational cycles of LST-1. The result for P2 proves the excellent detection capabilities of the LST-1 compared to other existing IACTs: we obtained a doubled significance when compared to MAGIC Collaboration (2020) in less observational time and with a single telescope.

The absence of emission from P1 above 15 GeV was already reported by Ackermann et al. (2013) and later by MAGIC Collaboration (2020). As a further comparison, we obtained the Fermi-LAT phaseograms using a circular extraction region of $3^{\circ}$ radius for different energy ranges, shown in Fig. 2, to see how the two peaks evolve. The significance of both peaks and the Bridge is above $100\sigma$ up to 10 GeV, but it decreases to $5.7\sigma$, $23\sigma$ and $4\sigma$ above 10 GeV for P1, P2, and the Bridge, respectively. Above 15 GeV, only P2 is detected with a significance of $10.7\sigma$, while P1 and the Bridge both show excesses at $1\sigma$ level. The results obtained with the LST-1 are in line with these findings. However, it is possible to observe some hint of excess (significance above $1.5\sigma$) from P1 and the Bridge regions in the LST-1 sample, as shown in Fig. 1. These features will be investigated in the future when additional LSTs are operative.

## 3.2. Morphology of P2

We studied the shape of the second peak in the LST-1 phaseogram, fitting both a symmetric Gaussian and a symmetric Lorentzian model in phase. We chose symmetric models since

no clear evidence of asymmetries is observed in P2 (Abdo et al. 2010). We added a constant representing the background to each profile so that the final models could be written as

$$f_{\text{Gauss}}(\phi) = A_G \exp\left(-\frac{(\phi - \mu)^2}{2\sigma^2}\right) + \text{const}_{\text{Gauss}} \quad (1)$$

for the Gaussian, and

$$f_{\text{Lor}}(\phi) = A_L \left[ 1 + \left( \frac{\phi - x_0}{2\gamma} \right)^2 \right]^{-1} + \text{const}_{\text{Lor}} \quad (2)$$

for the Lorentzian. In both cases, $\phi$ represents the value of the rotational phase.

We fitted these models in different reconstructed energy bins. We analysed the reconstructed energy range [15, 65] GeV, where the signal of P2 is concentrated. We further divided it into two logarithmically spaced bins, [15, 31] GeV and [31, 65] GeV, and repeated the analysis in both. We did not include data above 65 GeV due to the lack of signal from P2 at those energies.

The best-fit results for the full band and the two energy bins are reported in Table 1 for both the Gaussian and the Lorentzian profiles, along with their statistical errors. All the statistical uncertainties from now on in the paper will be reported at $1\sigma$ level. The full width at half maximum (FWHM) has been chosen as a measure of the peak width and it has been computed as $FWHM_{\text{Gauss}} = 2\sigma\sqrt{2\log 2}$ for the Gaussian and as $FWHM_{\text{Lor}} = 2\gamma$ for the Lorentzian. Henceforth, the notation ‘log’ refers to the natural logarithm.

We found, for both models, that the position of the second peak of Geminga does not significantly change in the two studied energy bins. The value is compatible with that reported in Ceribella (2021) within the statistical errors, but slightly differs from the best-fit values of Abdo et al. (2010). The same conclusions can be drawn for the peak width results.

The results obtained for the Gaussian and Lorentzian models are consistent with each other. To assess the goodness of fit of both models, we reported the associated p-value. We also computed the Akaike information criterion (Akaike 1974) defined as $AIC = 2k - 2\log(L)$, where k and L represent the number of model parameters and the maximum likelihood of the model. The AIC can be used to compare the results of two non-nested models and determine if one of the two is preferred. In general, lower coefficient values indicate a better agreement between the fit and the data. For the LST-1 sample, the Gaussian profile shows slightly lower values of the AIC. We computed the difference $\Delta(AIC)$ for both the full band and the energy bins, obtaining $\Delta(AIC) = (9.6, 2.2, 9.2)$ for the broad-band, the first and second bin, respectively.

As an additional cross-check, we fitted the phaseograms obtained with the Fermi-LAT sample described in Sect. 2.2 for energies above 15 GeV (considered as the ‘broad-band’ Fermi sample) and in the same energy bins as for the LST-1 analysis, i.e. [15, 31] GeV and [31, 65] GeV. The results obtained, shown in Table 1, are consistent with the LST-1 ones within the statistical uncertainties for both the Gaussian and the Lorentzian profiles. The fit in the second energy bin did not converge well, probably due to the low statistics of the Fermi-LAT sample above 30 GeV, and for this reason we did not report the results. Also in this case, there is no significant preference between Gaussian and Lorentzian profiles, since both the p-values and the AIC in both the broad-band and the [15, 31] GeV fits are similar.
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Phaseogram of Geminga obtained from the analysis of 16.6 years of Fermi-LAT data using a circular extraction region with a 3° radius for four different energy ranges, indicated on top of each plot, and using the same colour code for the phase regions as Fig. 1.</figcaption></figure>

<table><caption>**Table 1.** Best-fit results for the P2 mean position and width (FWHM) for the symmetric Gaussian and Lorentzian profiles.</caption><thead><tr><th rowspan="2">Instrument</th><th colspan="5">Gaussian</th><th colspan="3">Lorentzian</th></tr><tr><th>Energy (GeV)</th><th>μ</th><th>FWHM</th><th>p-value</th><th>AIC</th><th>x<sub>0</sub></th><th>FWHM</th><th>p-value</th><th>AIC</th></tr></thead><tbody><tr><td rowspan="3">CTAO LST-1</td><td>[15, 31]</td><td>0.595 ± 0.002</td><td>0.051 ± 0.004</td><td>0.91</td><td>51.3</td><td>0.597 ± 0.002</td><td>0.043 ± 0.005</td><td>0.62</td><td>60.9</td></tr><tr><td>[31, 65]</td><td>0.597 ± 0.004</td><td>0.05 ± 0.01</td><td>0.74</td><td>57.0</td><td>0.597 ± 0.004</td><td>0.04 ± 0.01</td><td>0.66</td><td>59.1</td></tr><tr><td>Full band</td><td>0.596 ± 0.002</td><td>0.050 ± 0.005</td><td>0.60</td><td>61.3</td><td>0.597 ± 0.002</td><td>0.041 ± 0.006</td><td>0.28</td><td>70.5</td></tr><tr><td>Fermi-LAT</td><td>[15, 31]</td><td>0.594 ± 0.001</td><td>0.053 ± 0.003</td><td>0.74</td><td>197.9</td><td>0.594 ± 0.001</td><td>0.050 ± 0.004</td><td>0.36</td><td>201.7</td></tr><tr><td></td><td>>15</td><td>0.595 ± 0.001</td><td>0.047 ± 0.003</td><td>0.67</td><td>224.5</td><td>0.595 ± 0.001</td><td>0.045 ± 0.008</td><td>0.69</td><td>222.7</td></tr><tr><td>MAGIC (Ceribella 2021)</td><td>[25, 100]</td><td>0.597 ± 0.005</td><td>0.061 ± 0.012</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

**Notes.** We show the results for both LST-1 and Fermi-LAT, in the two energy bins and the full band, with the associated statistical uncertainties and the correspondent AIC and p-value. As a comparison, the MAGIC results from Ceribella (2021) are also included in the table. The results for the [31, 65] GeV fit of the Fermi-LAT sample are not reported since the fit did not converge well.

### 3.3. Spectral energy distribution of P2

In order to obtain the SED of the second peak, we performed a forward folding fit using Gammapy in the energy range [20, 95] GeV, considering five bins per decade, using a power law (PL) model:

$$ \frac{dN}{dE} = f_0 \left( \frac{E}{E_0} \right)^{-\Gamma} \qquad (3) $$

where $f_0$ represents the flux normalisation, $\Gamma$ is the spectral index, and $E_0$ is the reference energy.

The choice for the lower edge of the fitting range is connected to the estimated energy threshold of the analysis. We chose the upper edge of our energy binning closest to 100 GeV

to assess the possibility of obtaining a spectral point at energies above 65 GeV. For more details, refer to Sect. 3.3.1. We set the reference energy to $E_0 = 14.3$ GeV, the decorrelation energy computed with Gammapy, to have the lowest correlation between the parameters as well as the lowest statistical error on the flux normalisation.

The fit results with the power law model on P2 are reported in Table 2, together with their statistical uncertainties. The spectrum observed by the LST-1 is slightly harder than that of MAGIC Collaboration (2020), but the two spectral indices are compatible if the systematic uncertainties are considered (see Sect. 4).

The flux points were calculated assuming the power law model fitted to the data, fixing the index and leaving free the nor-
<table><caption>**Table 2.** Best-fit results of the spectral fitting of P2 with a power law model.</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Γ</td><td>4.5 ± 0.4</td></tr><tr><td>f<sub>0</sub></td><td>(9.99 ± 0.75) · 10<sup>-8</sup> cm<sup>-2</sup> s<sup>-1</sup> TeV<sup>-1</sup></td></tr></tbody></table>

**Notes.** Each parameter is shown with its statistical uncertainty. The reference energy is fixed to $E_0 = 14.3$ GeV.

<table><caption>**Table 3.** Flux points of the P2 SED obtained assuming the best-fit PL model of Table 2.</caption><thead><tr><td>$E_{\text{ref}}$ (GeV)</td><td>$E_{\text{min}}$ (GeV)</td><td>$E_{\text{max}}$ (GeV)</td><td>Flux (TeV·cm<sup>-2</sup>·s<sup>-1</sup>)</td></tr></thead><tbody><tr><td>20</td><td>17</td><td>24</td><td>(9.2 ± 1.1) · 10<sup>-12</sup></td></tr><tr><td>28</td><td>24</td><td>33</td><td>(3.5 ± 0.4) · 10<sup>-12</sup></td></tr><tr><td>40</td><td>33</td><td>47</td><td>(1.6 ± 0.4) · 10<sup>-12</sup></td></tr><tr><td>56</td><td>47</td><td>67</td><td>(1.2 ± 0.4) · 10<sup>-12</sup></td></tr></tbody></table>

**Notes.** The upper limit on the last energy bin is not shown.

malisation in each bin, a procedure similar to the one explained in Sect. 2.2. Both the LST-1 best-fit model and the flux points are depicted in Fig. 4, and the latter are also reported in Table 3. The signal from P2 is detected up to ~65 GeV, in the last energy bin we computed an upper limit. This is compatible with what has been found when studying the phaseogram (see Sect. 3.1).

### 3.3.1. Energy threshold of the analysis and spectral binning

To evaluate the analysis threshold we fitted with a Landau profile the energy distribution of the MCs employed in the data reduction process, in order to extract the position of the peak energy. The MCs are simulated with PL spectrum with index $\Gamma = 2$, so we weighted the distribution for the spectral index of Geminga ($\Gamma = 4.5$). The energy threshold for the lowest-zenith MC determines the energy threshold of the analysis and it was found to be ~10 GeV in true energy. This value is lower than the one reported in Abe & Abe (2024) for the Crab pulsar study, due to the softness of the Geminga spectrum when compared to the Crab, and comparable with the performance of MAGIC's Sum-Trigger-II shown in MAGIC Collaboration (2020). However, due to Geminga's spectrum being very steep, the bias, i.e. the mean difference between the reconstructed and the true energy, is considerably high around the threshold. Indeed, when producing the one-dimensional histogram of the MC distribution as a function of the reconstructed energy instead of the true energy, we found the peak was located around 20 GeV. This implies the energy threshold of 10 GeV in true energy corresponds to ~20 GeV in reconstructed energy.

To better visualise the results, we produced the two-dimensional weighted histogram of the reconstructed energy versus the true energy, projected onto the true energy axis. This is the same as producing the one-dimensional weighted histogram of the MCs. The 2-D histogram, nevertheless, allows us to understand from which true energy the reconstructed events come from. We normalised the histogram to have a rate on the z-axis and plotted it for the two reconstructed energy bins used in the morphological study, i.e. [15, 31] and [31, 65] GeV. Figure 3 depicts the results. The Monte Carlo data chosen for Fig. 3 were simulated at $Zd = 10^{\circ}$, the lowest zenith for the MC produc-

tion, in order to obtain the plots for the lowest energy threshold possible.

For the energy binning of the SpectrumDataset Gammapy object representing the spectrum, we defined the true energy axis $E_{\text{true}}$ with 100 logarithmically spaced bins between 3 GeV and 50 TeV and the reconstructed energy axis $E_{\text{reco}}$ with 40 logarithmically spaced bins between 10 GeV and 10 TeV. The range for the spectral fit was defined in $E_{\text{reco}}$, due to the requirements of Gammapy, and it was chosen as a subset of the correspondent axis. We chose to start the subset from the edge closest to 20 GeV and end it at the edge closest to 100 GeV. This roughly corresponds to ~10–100 GeV in true energy. For the flux point calculation, we set the first edge of the energy binning at $E_{\text{reco}} \approx 17$ GeV to ensure the first flux point is around 20 GeV.

### 3.4. Joint Fermi-LAT and LST-1 SED of P2

To characterise the gamma-ray emission in a larger energy range we included in the fit the Fermi-LAT flux points computed with the sample described in Sect. 2.2. This is needed, in particular, to determine if a spectral cut-off exists and to verify if the emission at different energies is generated via the same mechanisms.

We performed a joint fit of LST-1 and Fermi-LAT data between 100 MeV and 100 GeV testing two spectral models. The first was a power law with an exponential cut-off (PLEC), and the second was a power law with a sub-exponential cut-off (PLSEC). The same mathematical law describes both models:

$$ \frac{dN}{dE} = f_0 \left(\frac{E}{E_0}\right)^{-\Gamma} \exp \left[ - \left(\frac{E}{E_c}\right)^{\beta} \right] = f_0 \left(\frac{E}{E_0}\right)^{-\Gamma} \exp \left[ - (\lambda E)^{\beta} \right], \quad (4) $$

where $\beta = 1$ in the case of PLEC, while $\beta < 1$ for the PLSEC. $E_0$ is the reference energy and $E_c$ is the cut-off energy. We also report the formulation of the equation used in Gammapy, which uses the reciprocal of the cut-off energy $\lambda = 1/E_c$.

The results of the joint fit are described in Table 4 and the plot is shown in Fig. 4, along with the MAGIC Collaboration (2020) points for comparison. The value of the $\beta$ coefficient for the PLSEC model is compatible with that obtained by MAGIC Collaboration (2020) considering the statistical uncertainties.

Since the PLSEC and PLEC are nested models, Wilks theorem (Wilks 1962) allows us to use of the likelihood ratio test to determine which is preferred. The PLEC model was found to be disfavoured with $-2\Delta \log \mathcal{L} = 323$. Considering a chi-square distribution with one degree of freedom, it corresponds to a rejection level of $TS \approx 18\sigma$. This is in agreement with the findings already presented in the second pulsar catalogue (2PC) of Fermi-LAT (Abdo et al. 2013) and later confirmed in the 3PC (Smith et al. 2023), which found a preference in pulsar spectra for a softer sub-exponential cut-off rather than the classical exponential. This is also compatible with what has been found by MAGIC Collaboration (2020). However, it is clear from the plot that the PLSEC still does not match well with the LST-1 flux points, even though it is preferred over the PLEC. To better investigate the preference for a non-curved model over a curved one, we performed a test to assess the presence of curvature (see Sect. 3.5).

### 3.5. Testing for curvature in the spectrum

We used the Fermi-LAT flux points to verify the presence of a possible curvature in the spectrum at energies above 10 GeV. To do so, we fitted the LST-1 and Fermi points using a power
<figure><img src="image_4.png" /><figcaption>Fig. 3. Two-dimensional weighted (considering Geminga's spectral index Γ = 4.5) histogram of the reconstructed energy versus the true energy, projected onto the true energy axis for the two reconstructed energy bins used for the morphological study of P2, [15, 31] GeV (left) and [31, 65] GeV (right). The Monte Carlo data used for the plots were produced at Zd = 10°. The dashed line represents the equivalence between the true energy E<sub>true</sub> and the reconstructed energy E<sub>reco</sub>. The z-axis is in units of rate, i.e. events per second.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4. Joint LST-1 (squares) and Fermi-LAT (circles) data samples of P2, along with the best-fit results of both the power law with an exponential cut-off (PLEC, dotted line) and the power law with sub-exponential cut-off (PLSEC, dashed line). The power law fit of the LST-1 only points (orange squares) is shown together with its statistical 1σ uncertainty band (solid line and shaded area) and the systematics uncertainty band (dash-dotted area), considering both the systematics on the index and the flux normalisation. The MAGIC Collaboration (2020) points are depicted as triangles for comparison. The horizontal error bars represent the width of the energy bins.</figcaption></figure>

law and a log parabola, the latter described by the following parametrisation:

$$ \frac{dN}{dE} = f_0 \left( \frac{E}{E_0} \right)^{-a-b \log(E/E_0)} \quad (5) $$

where the parameter b represents the curvature index, and it is equal to zero for a non-curved model (i.e. a simple power law). We fitted both models to the joint dataset imposing a minimum energy of the Fermi-LAT points of 10, 15 and 20 GeV. We also performed this test on the LST-1-only sample, using 20 and 25 GeV as minimum energies for the fit.

When considering the joint Fermi-LAT and LST-1 sample, a log parabola with concave (negative) curvature (i.e. positive b) is strongly favoured over the power law when setting a minimum energy of 10 GeV and 15 GeV, with ΔTS = 10σ for the former and ΔTS = 8.9σ for the latter. In both cases, the best-fit curvature parameter is b = 3.3 ± 0.6. The preference for a curved log parabola model decreases to a 4.6σ level for

$E_{min} = 20$ GeV, for which b = 1.95 ± 0.99 is derived. We also tested $E_{min} = 25$ GeV and obtained again a slight preference for the log parabola at 2.6σ level, with best-fit curvature parameter b = 0.9 ± 0.3, almost consistent with zero. The overall preference for a curved log parabola model is probably related to the larger number of Fermi-LAT points when compared to those of LST-1, especially in the range 10–20 GeV that is not covered by the latter instrument. This preference tends to decrease as $E_{min}$ increases.

For the Fermi-LAT only sample, instead, we obtained a 3.5σ preference for the log parabola and a best fit b = 1.0 ± 0.3 when setting $E_{min} = 10$ GeV, and only a 1σ preference for a log parabola with b always consistent with zero in the case of $E_{min} = 15, 20$ GeV. No curvature seems to be detected at higher energies and it is likely due to the low statistics of the sample in the correspondent energy range.

Similarly it happens in the case of the LST-1 only sample, for which the results obtained fitting a log parabola for both minimum energies return a value of b which is consistent with zero and a $-2\Delta\log L \ll 1$ when compared to the power law. This implies no statistical preference between the two models and no curvature is detected in the LST-1 sample. However, a non-detection does not necessarily mean the curvature is absent for a wider energy range. With the current LST-1 sample, we cannot discard either possibility, but it may be possible to understand more about it by collecting more data to lower the statistical uncertainty and extending it to higher energies.

## 4. Evaluation of the systematic effects

The best-fit result for the PL fitting of the P2 spectrum we found was slightly harder than the spectrum reported by MAGIC Collaboration (2020). To assess the compatibility of MAGIC's results with ours, we examined the contribution of the systematic effects in the analysis of the LST-1-only sample. Different factors can contribute to the systematic uncertainties, with a more or less significant impact on the final result, and we investigated the most relevant ones. This is the first time such a detailed study of the systematics is performed for LST-1 analysis.
<table><caption>**Table 4.** Best-fit results for the joint fit of Fermi-LAT and LST-1 data, for both tested models.</caption><thead><tr><th>Model</th><th>f<sub>0</sub> (TeV<sup>-1</sup> cm<sup>-2</sup> s<sup>-1</sup>)</th><th>E<sub>0</sub> (GeV)</th><th>Γ</th><th>λ (GeV<sup>-1</sup>)</th><th>E<sub>c</sub> (GeV)</th><th>β</th><th>-2 log L</th></tr></thead><tbody><tr><td>PLEC</td><td>(3.610 ± 0.002) · 10<sup>-4</sup></td><td>1.00</td><td>1.079 ± 0.005</td><td>0.346 ± 0.003</td><td>2.890 ± 0.003</td><td>1</td><td>430.3</td></tr><tr><td>PLSEC</td><td>(6.1 ± 0.3) · 10<sup>-4</sup></td><td>1.00</td><td>0.85 ± 0.02</td><td>0.76 ± 0.05</td><td>1.32 ± 0.05</td><td>0.73 ± 0.02</td><td>107.3</td></tr></tbody></table>

**Notes.** The critical energy $E_c$ is computed from its reciprocal $\lambda$ estimated with Gammapy.

## 4.1. Different MC efficiency

The first factor we studied was the setting of different energy-dependent selection efficiencies of the MC for both alpha and gammaness. The value of 70% is usually chosen as the standard because it is a good compromise between having a large number of events classified as gamma-originated and a reliable gamma-hadron separation. To test whether and how the choice of efficiency impacts the final results, we computed the spectral parameters in a grid of alpha efficiencies of 70% and 90% and gammaness efficiencies of 40%, 70%, and 90%.

We found that both the spectral index and the flux normalisation change within 5% of the reference value, with a maximum positive absolute variation of $\Gamma$ of +0.13 and a maximum negative absolute variation of -0.02. We can conclude that this effect is subdominant to the statistical uncertainties.

## 4.2. Mismatches in the reconstruction

Another important systematic effect that must be taken into account is related to the mismatches between the MC and the observed data, which lead to a wrong reconstruction of the events. One of the most common ones is the energy scale mismatch considered here. Non-ideal conditions during the observations, such as dirty mirrors or lower atmospheric transparency, affect the total amount of recorded light by the system, introducing a bias in the reconstruction, in particular, of the energy. The proper way to simulate this effect would be to have dedicated MC productions in which an energy bias is introduced. Still, this method is time-consuming and requires significant resources. Defining the energy scaling as $\epsilon$, we can say that, for small values of $\epsilon$, the effect of the energy bias can be approximately reproduced by shifting the migration matrix along the true energy axis.

We considered $\epsilon = \pm 15\%$, which is the typical value used for the studies of energy scale systematics for IACTs (Aleksic et al. 2016). We produced the scaled IRFs in the two cases and derived the spectrum again. We found that the spectral index is not significantly affected by the energy scaling, while the flux normalisation changes even up to ~55% as a consequence of the shift of the spectrum in the energy axis. In detail, the absolute positive variation of $\Gamma$ is +0.04 and the negative is -0.08. The overall effect is a rigid shift of the PL to the right or the left depending on the scaling.

The same type of systematic uncertainty has been considered in MAGIC Collaboration (2020), for which they reported a ~1% contribution for the spectral index, corresponding to an absolute variation of ~0.06. If we only consider this effect, the results of LST-1 are still not consistent with MAGIC's.

## 4.3. Different Zd cut of the sample

The choice of the zenith angle cut of the data impacts the overall energies of the events: the larger the applied zenith cut, the higher the energy threshold of the data sample will be. We performed the analysis using a zenith cut of $Zd < 50^{\circ}$ to have the largest possible sample, but ~90% of the observation time (~54 hours out of 60) data were taken at $Zd < 25^{\circ}$, assuring the lowest possible energy threshold for the sample. In Table 5, the observational hours for each of the Zd cuts are shown.

<table><caption>**Table 5.** Total observation time of the LST-1 sample for the different maximum zenith distance (Zd) cuts applied for the analysis of the systematic effects.</caption><thead><tr><th>Max Zd (deg)</th><th>Total time (h)</th><th>Max Zd (deg)</th><th>Total time (h)</th></tr></thead><tbody><tr><td>50</td><td>60.1</td><td>25</td><td>54.1</td></tr><tr><td>40</td><td>59.5</td><td>20</td><td>42.8</td></tr><tr><td>35</td><td>58.7</td><td>15</td><td>26.8</td></tr><tr><td>30</td><td>56.6</td><td></td><td></td></tr></tbody></table>

We obtained the spectrum for different values of the zenith cut of the sample and estimated how much the choice of Zd influences $\Gamma$ and $f_0$. The most significant variation is obtained when setting $Zd < 20^{\circ}$. In this case, we obtain an absolute negative variation in the index of -0.5 and a relative positive variation in the normalisation of +12%. The other variations obtained with this test are much smaller.

## 4.4. Different intensity cut of the sample

The last investigated systematic effect is related to the choice of the intensity cut applied to both simulated and observed data. This cut is usually applied to guarantee an overall good match between simulations and observations, and the standard value adopted for LST-1 analysis is an image intensity above 50 p.e. as shown in CTAO-LST Project (2023). We tried a series of different cut values, from a minimum intensity of 20 p.e. to a maximum of 70 p.e., and computed the spectral results in each case.

In this case, the positive and negative absolute variations of the spectral index are +0.2 and -0.3, respectively. The uncertainties are around 20% for the normalisation.

## 4.5. Total contribution

To compute the total contribution of the systematic uncertainties on the spectral index, we assumed the correlation between all the examined effects to be low. In this way, we could sum in quadrature all the maximum positive and negative variations to obtain the total positive and negative uncertainty. We derived a final result for the spectral index of

$$ \Gamma = (4.5 \pm 0.4_{\text{stat}})^{+0.2_{\text{sys}}-0.6_{\text{sys}}} $$

This result is compatible with the one reported by MAGIC Collaboration (2020), considering the statistical uncertainty of the MAGIC result. This is the first time another facility has cross-checked the MAGIC Collaboration's result and it confirms again the potential of the LSTs.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Adaptation of Fig. 10 of Harding et al. (2021) showing the two models derived in the paper for the SC emission of Geminga, together with the LST-1 P2 SED obtained in this work (same as Fig. 4). The Fermi-LAT and MAGIC points from the analysis presented in this work and MAGIC Collaboration (2020), respectively, are also reported for comparison. Note that the Fermi-LAT points represent the phase-averaged flux, i.e. they take into account the emission from both P1 and P2.</figcaption></figure>

The systematic uncertainties on the normalisation are much larger, $\Delta f_{0,\text{sys}}^+ = +6 \cdot 10^{-8} \text{ cm}^{-2} \text{ s}^{-1} \text{ TeV}^{-1}$ and $\Delta f_{0,\text{sys}}^- = -5 \cdot 10^{-8} \text{ cm}^{-2} \text{ s}^{-1} \text{ TeV}^{-1}$. This is expected since the spectrum of Geminga is very steep and small fluctuations of the spectral index can induce large variations of the normalisation. The total systematic uncertainty band, considering both the uncertainty on the index and normalisation, is depicted in Fig. 4 and in Fig. 5 as the dash-dotted area.

## 5. Discussion and conclusions

In this paper, we report the detection of the Geminga pulsar by the LST-1 of CTAO. The signal from the second peak of the phaseogram, P2, is detected with a significance of $12\sigma$, achieved with 60 hours of good-quality data. The first peak, P1, remains undetected, even though a $2.6\sigma$ excess has been observed. We also observed an excess of $1.8\sigma$ for the first inter-peak (Bridge) region.

The detection of P1 with the LST-1 alone is still challenging. Considering the evolution of the significance with time we derived in Sect. 3.1, more than 200 hours would be needed to achieve a $5\sigma$ detection. However, with the deployment of the full array of four LSTs at CTAO-North, it will be possible to reduce this time and detect P1, and possibly the Bridge, in the near future. Using the latest ProD5 IRFs of CTAO (Cherenkov Telescope Array Observatory & Cherenkov Telescope Array Consortium 2021) and an extrapolated power law spectrum from the Fermi-LAT data above 5 GeV, we estimated around 30 hours will be needed for a $5\sigma$ detection of the first peak with the full array of LSTs.

We performed a morphological study of the second peak by testing two models, the symmetric Gaussian and the symmetric Lorentzian. We found compatible peak mean positions and FWHMs between the two profiles and no energy-dependent evolution in the two logarithmically spaced energy bins considered, [15, 31] GeV and [31, 65] GeV. The Gaussian model, in general, returns better goodness-of-fit values, highlighting a slight preference for this profile. More tests on the pulse morphology will be carried out after new observations, possibly extending the study to P1 and the Bridge, to further investigate the geometry of the peaks.

We then extracted the SED of P2 in the [20, 95] GeV range, obtaining four flux points and an upper limit in the last energy bin. As also observed with the study of the phaseogram, the signal from P2 is reconstructed up to 65 GeV. The best-fit result yielded a power law with spectral index $\Gamma = (4.5 \pm 0.4)_{\text{stat}}^{+0.2_{-0.6}}_{-0.2_{+0.6}}$, compatible with the previous result of MAGIC Collaboration (2020), and a normalisation at 14.3 GeV of $f_0 = (9.99 \pm 0.75)_{\text{stat}}^{+6_{-5}}_{-5} \cdot 10^{-8} \text{ cm}^{-2} \text{ s}^{-1} \text{ TeV}^{-1}$. For the first time, a result on a pulsar other than the Crab is being cross-checked and confirmed by two different IACTs. This is important in the field of pulsars detected by IACTs, since only a few sources have been detected and it is not straightforward to perform cross-checks between different instruments.

We also analysed 16.6 years of Fermi-LAT data to perform a joint fit with the LST-1 sample and extend the Geminga spectrum down to 100 MeV. Two models, a power law with exponential or sub-exponential cut-off, were tested and we found that the exponential cut-off model can be rejected at an $18\sigma$ level, in line with previous findings by Fermi-LAT and MAGIC. However, even if the power law with sub-exponential cut-off returns a better fit, it does not correctly describe the spectrum obtained with LST-1.

For this reason, using the Fermi flux points, we tested for the presence of curvature in the high-energy end of the spectrum. We performed a likelihood ratio test between a power law and a log parabola spectral model using different minimum energies for the Fermi-LAT data, and we repeated the test on the LST-1-only sample. For the latter, no curvature is detected, probably due to the limited statistics and sampled energy range. When considering the joint dataset, the log parabola is found to be statistically favoured over the power law. However, the preference tends to decrease as the minimum energy set as the starting point of the fit increases. Nevertheless, the fit results may be affected by the absence of LST-1 flux points below 20 GeV or by the low statistics of both samples at the highest energies.

After the early measurements by Fermi-LAT, theories based on the outer gap (OG) concept were favoured (Cheng et al. 1986; Romani 1996). To explain the gamma-ray emission from pulsars, these theories assumed that particles were accelerated to relativistic velocities inside charge-depleted regions high up in the magnetosphere. However, in addition to the uncertainties and approximations used by these models (see Viganò et al. (2015a) for a discussion), particle-in-cell (PIC) and force-free electrodynamics simulations signalled that instead of gaps, the most likely location of the acceleration of the gamma-ray-emitting particles is the so-called Y-point, a localised region at or beyond the light cylinder, see Cerutti (2019) and references therein for details. Also, particularly after the second Fermi-LAT pulsar catalogue (Abdo et al. 2013), it was clear that the variety of the Fermi-LAT spectra, even without considering VHE tails, was already beyond what curvature-only models can describe. Viganò et al. (2015b) proposed a synchro-curvature (SC) approach to tackle this variety of spectra. This was used to provide fits of the spectra of gamma and X-ray pulsars (Torres 2018), and when coupled to light curves predictions, to describe their global properties as well (Íñiguez-Pascual et al. 2024). Harding et al. (2021), Kalapotharakos et al. (2018), Brambilla et al. (2018) and Cerutti et al. (2025), among a few others, introduced PIC-motivated models. In all these cases, though, the broad-band emission, particularly beyond tera-electronvolt energies, seems to require Compton components.

Harding et al. (2021) proposed to describe the optical to X-ray emission with synchrotron radiation, SC (as used above) up to ~100 GeV and inverse Compton scattering (ICS) above
tera-electronvolt energies. They also include a synchrotron-self-Compton (SSC) component, but, except for the Crab pulsar, the SC emission completely outshines it. The ICS component at tera-electronvolt energies, instead, could in principle be detected. The Harding et al. (2021) model comprises several parameters, such as the inclination angle $\alpha$ between the rotational axis and the magnetic field, the viewing angle $\zeta$ between the pulsar and the observer, and the low and high values of the accelerated electric field, the latter expressed through the parameters $R_{\text{acc}}^{\text{low}} = eE_{\parallel}^{\text{low}}/\text{mc}^2$ and $R_{\text{acc}}^{\text{high}} = eE_{\parallel}^{\text{high}}/\text{mc}^2$. The values of these parameters are adjusted to match the observations, in particular $E_{\parallel}^{\text{low}}$ and $E_{\parallel}^{\text{high}}$ are chosen to reproduce the hard X-rays to gigaelectronvolt data and $\zeta$ is fixed to the value that best matches the Fermi-LAT and VHE data. After setting the values of those parameters, the ICS and the total SED are computed.

For the specific case of Geminga, the SED model was obtained for $\alpha = 75^{\circ}$, $\zeta = 55^{\circ}$ and for two different sets of values for the low and high accelerating electric field, the first one with $R_{\text{acc}}^{\text{low}} = 0.04$ and $R_{\text{acc}}^{\text{high}} = 0.15$ and the second with $R_{\text{acc}}^{\text{low}} = R_{\text{acc}}^{\text{high}} = 0.15$. As one can see from Fig. 10 of Harding et al. (2021), the dominant component of the high-energy emission of Geminga is that of SC radiation. Moreover, the predicted ICS is too faint to be detected, even considering CTAO-North's sensitivity. This is caused by the larger radius of the light cylinder due to the lower rotational speed of Geminga when compared to the younger pulsars (i.e. Vela, Crab) considered in Harding et al. (2021), which implies a larger distance between the current sheet and the neutron star producing the X-ray photon seeds for the ICS. In this context, the spectrum observed both by MAGIC and LST-1 is a continuous tail-like extension of the Fermi-LAT one and is interpreted as curvature radiation emission, without the need for additional components. We depicted in Fig. 5 the LST-1 SED together with the SC component of Geminga predicted by this model for the two cases of $R_{\text{acc}}^{\text{low}}$, $R_{\text{acc}}^{\text{high}}$ considered in the reference, as explained above. The model describes the emission of both peaks together, so we included in the same Figure the Fermi-LAT flux points of the phase-averaged analysis. In the case of MAGIC and LST-1, instead, we considered the flux points of P2 only, since the emission above 15 GeV is fully dominated by the second peak. The Harding et al. (2021) model also predicts the absence of signal from P1 above 20 GeV, which is compatible with what has been found in this work. All taken into account, this model could explain some phenomenology of the Geminga emission from the Fermi-LAT to the LST-1 energy range, even though the match with the LST-1 SED and the higher energy flux points of Fermi-LAT could be improved. With the present sample, it is not possible to completely rule out this or other theoretical interpretations and many open questions remain. Further studies, possibly with more LSTs, are necessary to validate or reject the model.

Overall, we can conclude that the results obtained on the Geminga pulsar prove the LST-1 to be an excellent telescope for pulsar observations at the upper end of their spectra, also considering its partial overlap with the Fermi-LAT energy range. With this study, we were able to prove that the threshold of the telescope for sources with very steep spectra can be as low as 10 GeV, even though a large bias is observed at energies below 20 GeV. This issue will be reduced when at least another LST becomes operational because the use of a stereoscopic system improves the reconstruction of the events. With the full array of LSTs, the sensitivity will be improved and more tests could be performed, such as investigating the behaviour of P1 and of the bridge as a function of energy, as well as a better determination of the P2 spectrum to verify the validity of the theoretical models.

**Acknowledgements.** We gratefully acknowledge financial support from the following agencies and organisations: Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), Fundação de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ), Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Fundação de Apoio à Ciência, Tecnologia e Inovação do Paraná - Fundação Araucária, Ministry of Science, Technology, Innovations and Communications (MCTIC), Brasil; Ministry of Education and Science, National RI Roadmap Project DO1-153/28.08.2018, Bulgaria; Croatian Science Foundation (HrZZ) Project IP-2022-10-4595, Rudjer Boskovic Institute, University of Osijek, University of Rijeka, University of Split, Faculty of Electrical Engineering, Mechanical Engineering and Naval Architecture, University of Zagreb, Faculty of Electrical Engineering and Computing, Croatia; Ministry of Education, Youth and Sports, MEYS LM2023047, EU/MEYS CZ.02.1.01/0.0/0.0/16_013/0001403, CZ.02.1.01/0.0/0.0/18_046/0016007, CZ.02.1.01/0.0/0.0/16_019/0000754, CZ.02.01.01/00/22_008/0004632 and CZ.02.01.01/00/23_015/0008197 Czech Republic; CNRS-IN2P3, the French Programme d'investissements d'avenir and the Enigmass Labex. This work has been done thanks to the facilities offered by the Univ. Savoie Mont Blanc - CNRS/IN2P3 MUST computing center, France; Max Planck Society, German Bundesministerium für Bildung und Forschung (Verbundforschung / ErUM), Deutsche Forschungsgemeinschaft (SFBs 876 and 1491), Germany; Istituto Nazionale di Astrofisica (INAF), Istituto Nazionale di Fisica Nucleare (INFN), Italian Ministry for University and Research (MUR), and the financial support from the European Union - Next Generation EU under the project IR0000012 - CTA+ (CUP C53C22000430006), announcement N.3264 on 28/12/2021: "Rafforzamento e creazione di IR nell'ambito del Piano Nazionale di Ripresa e Resilienza (PNRR)"; ICRR, University of Tokyo, JSPS, MEXT, Japan; JST SPRING - JPMJSP2108; Narodowe Centrum Nauki, grant number 2019/34/E/ST9/00224, Poland; The Spanish groups acknowledge the Spanish Ministry of Science and Innovation and the Spanish Research State Agency (AEI) through the government budget lines PGE2022/28.06.000X.711.04, 28.06.000X.411.01 and 28.06.000X.711.04 of PGE 2023, 2024 and 2025, and grants PID2019-104114RB-C31, PID2019-107847RB-C44, PID2019-104114RB-C32, PID2019-105510GB-C31, PID2019-104114RB-C33, PID2019-107847RB-C43, PID2019-107847RB-C42, PID2019-107988GB-C22, PID2021-124581OB-I00, PID2021-125331NB-I00, PID2022-136828NB-C41, PID2022-137810NB-C22, PID2022-138172NB-C41, PID2022-138172NB-C42, PID2022-138172NB-C43, PID2022-139117NB-C41, PID2022-139117NB-C42, PID2022-139117NB-C43, PID2022-139117NB-C44, PID2022-136828NB-C42, PDC2023-145839-100 funded by the Spanish MCIN/AEI/10.13039/501100011033 and "and by ERDF/EU and NextGenerationEU PRTR; the "Centro de Excelencia Severo Ochoa" program through grants no. CEX2019-000920-S, CEX2020-001007-S, CEX2021-001131-S; the "Unidad de Excelencia María de Maeztu" program through grants no. CEX2019-000918-M, CEX2020-001058-M; the "Ramón y Cajal" program through grants RYC2021-032991-I funded by MICIN/AEI/10.13039/501100011033 and the European Union "NextGenerationEU"/PRTR and RYC2020-028639-I; the "Juan de la Cierva-Incorporación" program through grant no. IJCC2019-040315-I and "Juan de la Cierva-formación" through grant JDC2022-049705-I. They also acknowledge the "Atracción de Talento" program of Comunidad de Madrid through grant no. 2019-T2/TIC-12900; the project "Tecnologías avanzadas para la exploración del universo y sus componentes" (PR47/21 TAU), funded by Comunidad de Madrid, by the Recovery, Transformation and Resilience Plan from the Spanish State, and by NextGenerationEU from the European Union through the Recovery and Resilience Facility; "MAD4SPACE: Desarrollo de tecnologías habilitadoras para estudios del espacio en la Comunidad de Madrid" (TEC-2024/TEC-182) project funded by Comunidad de Madrid; the La Caixa Banking Foundation, grant no. LCF/BQ/PI21/11830030; Junta de Andalucía under Plan Complementario de I+D+I (Ref. AST22_0001) and Plan Andaluz de Investigación, Desarrollo e Innovación as research group FQM-322; Project ref. AST22_00001_9 with funding from NextGenerationEU funds; the "Ministerio de Ciencia, Innovación y Universidades" and its "Plan de Recuperación, Transformación y Resiliencia"; "Consejería de Universidad, Investigación e Innovación" of the regional government of Andalucía and "Consejo Superior de Investigaciones Científicas", Grant CNS2023-144504 funded by MICIU/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR, the European Union's Recovery and Resilience Facility-Next Generation, in the framework of the General Invitation of the Spanish Government's public business entity Red.es to participate in talent attraction and retention programmes within Investment 4 of Component 19 of the Recovery, Transformation and Resilience Plan; Junta de Andalucía under Plan Complementario de I+D+I (Ref. AST22_00001), Plan Andaluz de Investigación, Desarrollo e Innovación (Ref. FQM-322). "Programa Operativo de Crecimiento Inteligente" FEDER 2014-2020 (Ref. ESFRI-2017-IAC-12), Ministerio de Ciencia y Innovación, 15% co-financed by Consejería de Economía, Industria, Comercio y Conocimiento del Gobierno de Canarias; the "CERCA" program and the grants 2021SGR00426 and 2021SGR00679, all funded by the Generalitat de Catalunya; and the European Union's NextGenera-
tionEU (PRTR-C17.11). This research used the computing and storage resources provided by the Port d'Informació Científica (PIC) data center. State Secretariat for Education, Research and Innovation (SERI) and Swiss National Science Foundation (SNSF), Switzerland; The research leading to these results has received funding from the European Union's Seventh Framework Programme (FP7/2007-2013) under grant agreements No 262053 and No 317446; This project is receiving funding from the European Union's Horizon 2020 research and innovation programs under agreement No 676134; ESCAPE - The European Science Cluster of Astronomy & Particle Physics ESFRI Research Infrastructures has received funding from the European Union's Horizon 2020 research and innovation programme under Grant Agreement no. 824064.

**Author contribution:** G. Brunelli: project coordination, LST-1 pulsar analysis (phaseogram and SED), Fermi-LAT phaseograms, joint fit of LST-1 and Fermi-LAT, systematic error estimation, results discussion. A. Mas-Aguilar: analysis coordination, cross-check on LST-1 pulsar analysis (phaseogram and SED), cross-check on joint fit of LST-1 and Fermi-LAT, systematic error estimation, results discussion. G. Ceribella: Fermi-LAT analysis, ephemeris production, results discussion. P. K. H. Yeung: cross-check on LST-1 pulsar analysis (phaseogram and SED), cross-check of the Fermi-LAT phaseogram, results discussion. R. Lopez-Coto: PI of the observations, project and analysis coordination, results discussion. M. Lopez-Moya: Fermi-LAT analysis, results discussion. All corresponding authors have participated in the paper drafting and edition. The rest of the authors have contributed in one or several of the following ways: design, construction, maintenance and operation of the instrument(s) used to acquire the data; preparation and/or evaluation of the observation proposals; data acquisition, processing, calibration and/or reduction; production of analysis tools and/or related Monte Carlo simulations; discussion and approval of the contents of the draft.

## References

Abdo, A. A., Ackermann, M., Ajello, M., et al. 2010, ApJ, 720, 272

Abdo, A. A., Ajello, M., Allafort, A., et al. 2013, ApJS, 208, 17

Abdollahi, S., Acero, F., Ackermann, M., et al. 2020, ApJS, 247, 33

Abe, K., Abe, S., Aguasca-Cabot, A., et al. 2023, in Proceedings of 38th International Cosmic Ray Conference - PoS(ICRC2023), 444, 731

Ackermann, M., Ajello, M., Allafort, A., et al. 2013, ApJS, 209, 34

Aguasca-Cabot, A., Donath, A., Feijen, K., et al. 2023, https://doi.org/10.5281/zenodo.8033275

Akaike, H. 1974, IEEE Trans. Autom. Control, 19, 716

Aleksic, J., Alvarez, E. A., Antonelli, L. A., et al. 2012, A&A, 540, A69

Aleksic, J., Ansoldi, S., Antonelli, L., et al. 2016, Astropart. Phys., 72, 76

Aliu, E., Anderhub, H., Antonelli, L. A., et al. 2008, Science, 322, 1221

Ansoldi, S., Antonelli, L. A., Antoranz, P., et al. 2016, A&A, 585, A133

Bertsch, D. L., Brazier, K. T. S., Fichtel, C. E., et al. 1992, Nature, 357, 306

Brambilla, G., Kalapotharakos, C., Timokhin, A. N., Harding, A. K., & Kazanas, D. 2018, ApJ, 858, 81

Ceribella, G. 2021, Ph.D. Thesis, Technische Universität München

Cerutti, B. 2019, Rendiconti Lincei. Scienze Fisiche e Naturali, 30, 89

Cerutti, B., Figueiredo, E., & Dubus, G. 2025, A&A, 695, A93

Cheng, K. S., Ho, C., & Ruderman, M. 1986, ApJ, 300, 500

Cherenkov Telescope Array Observatory & Cherenkov Telescope Array Consortium 2021, https://doi.org/10.5281/zenodo.5499840

CTAO-LST Project, Abe, H., Abe, K., et al. 2023, ApJ, 956, 80

CTAO-LST Project, Abe, K., Abe, S., et al. 2024, A&A, 690, A167

Donath, A., Terrier, R., Remy, Q., et al. 2023, A&A, 678, A157

Ershov, A. A. 2021, Bull. Lebedev Phys. Inst., 48, 114

Faherty, J., Walter, F. M., & Anderson, J. 2007, Ap&SS, 308, 225

Fichtel, C. E., Hartman, R. C., Kniffen, D. A., et al. 1975, ApJ, 198, 163

Fomin, V. P., Stepanian, A. A., Lamb, R. C., et al. 1994, Astropart. Phys., 2, 137

Garcia, E., Vuillaume, T., & Nickel, L. 2022, arXiv e-prints [arXiv:2212.00120]

H.E.S.S. Collaboration (Aharonian, F., et al.) 2023, Nat. Astron., 7, 1341

Halpern, J. P., & Holt, S. S. 1992, Nature, 357, 222

Harding, A. K., Venter, C., & Kalapotharakos, C. 2021, ApJ, 923, 194

Íñiguez-Pascual, D., Torres, D. F., & Viganò, D. 2024, MNRAS, 530, 1550

Kalapotharakos, C., Brambilla, G., Timokhin, A., Harding, A. K., & Kazanas, D. 2018, ApJ, 857, 44

Li, T. P., & Ma, Y. Q. 1983, ApJ, 272, 317

López-Coto, R., et al. 2022, ASP Conf. Ser., 532, 357

López-Coto, R., Vuillaume, T., Moralejo, A., et al. 2024, https://doi.org/10.5281/zenodo.10604579

Luo, J., Ransom, S., Demorest, P., et al. 2019, PINT: High-precision pulsar timing analysis package, Astrophysics Source Code Library [record ascl:1902.007]

Luo, J., Ransom, S., Demorest, P., et al. 2021, ApJ, 911, 45

MAGIC Collaboration (Acciari, V.A., et al.) 2020, A&A, 643, L14

Mas, A., & Morcuende, D. 2024, https://doi.org/10.5281/zenodo.13385378

Romani, R. W. 1996, ApJ, 470, 469

Shearer, A., & Golden, A. 2002, in Neutron Stars, Pulsars, and Supernova Remnants, eds. W. Becker, H. Lesch, & J. Trümper, 44

Smith, D. A., Abdollahi, S., Ajello, M., et al. 2023, ApJ, 958, 191

Spir-Jacob, M., Djannati-Ataï, A., Mohrmann, L., et al. 2019, in Proceedings of 36th International Cosmic Ray Conference - PoS(ICRC2019), 358, 799

Torres, D. F. 2018, Nat. Astron., 2, 247

VERITAS Collaboration (Aliu, E., et al.) 2011, Science, 334, 69

Viganò, D., Torres, D. F., Hirotani, K., & Pessah, M. E. 2015a, MNRAS, 447, 2631

Viganò, D., Torres, D. F., Hirotani, K., & Pessah, M. E. 2015b, MNRAS, 447, 1164

Vuillaume, T., Garcia, E., & Nickel, L. 2024, https://doi.org/10.5281/zenodo.10649382

Wilks, S. S. 1962, Mathematical Statistics (John Wiley & Sons)

Wood, M., Caputo, R., Charles, E., et al. 2017, in 35th International Cosmic Ray Conference (ICRC2017), Int. Cosmic Ray Conf., 301, 824

Zanin, R., Abdalla, H., Abe, H., et al. 2021, in Proceedings of 37th International Cosmic Ray Conference - PoS(ICRC2021), 395, 005

1 Department of Physics, Tokai University, 4-1-1, Kita-Kaname, Hiratsuka, Kanagawa 259-1292, Japan

2 Institute for Cosmic Ray Research, University of Tokyo, 5-1-5, Kashiwa-no-ha, Kashiwa, Chiba 277-8582, Japan

3 INFN and Università degli Studi di Siena, Dipartimento di Scienze Fisiche, della Terra e dell'Ambiente (DSFTA), Sezione di Fisica, Via Roma 56, 53100 Siena, Italy

4 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, AIM, F-91191 Gif-sur-Yvette Cedex, France

5 FSLAC IRL 2009, CNRS/IAC, La Laguna, Tenerife, Spain

6 Departament de Física Quàntica i Astrofísica, Institut de Ciències del Cosmos, Universitat de Barcelona, IEEC-UB, Martí i Franquès, 1, 08028 Barcelona, Spain

7 Instituto de Astrofísica de Andalucía-CSIC, Glorieta de la Astronomía s/n, 18008 Granada, Spain

8 Department of Astronomy, University of Geneva, Chemin d'Ecogia 16, CH-1290 Versoix, Switzerland

9 INFN Sezione di Napoli, Via Cintia, ed. G, 80126 Napoli, Italy

10 INAF - Osservatorio Astronomico di Roma, Via di Frascati 33, 00040 Monteporzio Catone, Italy

11 Max-Planck-Institut für Physik, Boltzmannstraße 8, 85748 Garching bei München, Germany

12 INFN Sezione di Padova and Università degli Studi di Padova, Via Marzolo 8, 35131 Padova, Italy

13 Instituto de Astrofísica de Canarias and Departamento de Astrofísica, Universidad de La Laguna, C. Vía Láctea, s/n, 38205 La Laguna, Santa Cruz de Tenerife, Spain

14 Univ. Savoie Mont Blanc, CNRS, Laboratoire d'Annecy de Physique des Particules - IN2P3, 74000 Annecy, France

15 Universität Hamburg, Institut für Experimentalphysik, Luruper Chaussee 149, 22761 Hamburg, Germany

16 Graduate School of Science, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japan

17 IPARCOS-UCM, Instituto de Física de Partículas y del Cosmos, and EMFTEL Department, Universidad Complutense de Madrid, Plaza de Ciencias, 1. Ciudad Universitaria, 28040 Madrid, Spain

18 Faculty of Science and Technology, Universidad del Azuay, Cuenca, Ecuador

19 Centro Brasileiro de Pesquisas Físicas, Rua Xavier Sigaud 150, RJ 22290-180 Rio de Janeiro, Brazil

20 CIEMAT, Avda. Complutense 40, 28040 Madrid, Spain

21 University of Geneva - Département de physique nucléaire et corpusculaire, 24 Quai Ernest Ansernet, 1211 Genève 4, Switzerland

22 INFN Sezione di Bari and Politecnico di Bari, via Orabona 4, 70124 Bari, Italy

23 Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra, (Barcelona), Spain
24 INAF – Osservatorio Astronomico di Brera, Via Brera 28, 20121 Milano, Italy

25 Faculty of Physics and Applied Informatics, University of Lodz, ul. Pomorska 149-153, 90-236 Lodz, Poland

26 INAF – Osservatorio di Astrofisica e Scienza dello spazio di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy

27 Dipartimento di Fisica e Astronomia (DIFA) Augusto Righi, Università di Bologna, via Gobetti 93/2, I-40129 Bologna, Italy

28 Lamarr Institute for Machine Learning and Artificial Intelligence, 44227 Dortmund, Germany

29 INFN Sezione di Trieste and Università degli studi di Udine, via delle scienze 206, 33100 Udine, Italy

30 INAF – Istituto di Astrofisica e Planetologia Spaziali (IAPS), Via del Fosso del Cavaliere 100, 00133 Roma, Italy

31 Aix Marseille Univ, CNRS/IN2P3, CPPM, Marseille, France

32 University of Alcalá UAH, Departamento de Physics and Mathematics, Pza. San Diego, 28801 Alcalá de Henares, Madrid, Spain

33 INFN Sezione di Bari and Università di Bari, via Orabona 4, 70126 Bari, Italy

34 INFN Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy

35 Dipartimento di Fisica – Università degli Studi di Torino, Via Pietro Giuria 1, 10125 Torino, Italy

36 Palacky University Olomouc, Faculty of Science, 17. listopadu 1192/12, 771 46 Olomouc, Czech Republic

37 Dipartimento di Fisica e Chimica 'E. Segrè' Università degli Studi di Palermo, via delle Scienze, 90128 Palermo, Italy

38 INFN Sezione di Catania, Via S. Sofia 64, 95123 Catania, Italy

39 IRFU, CEA, Université Paris-Saclay, Bât 141, 91191 Gif-sur-Yvette, France

40 Port d'Informació Científica, Edifici D, Carrer de l'Albareda, 08193 Bellaterra (Cerdanyola del Vallès), Spain

41 INFN Sezione di Bari, via Orabona 4, 70125 Bari, Italy

42 Department of Physics, TU Dortmund University, Otto-Hahn-Str. 4, 44227 Dortmund, Germany

43 University of Rijeka, Department of Physics, Radmile Matejcic 2, 51000 Rijeka, Croatia

44 Institute for Theoretical Physics and Astrophysics, Universität Würzburg, Campus Hubland Nord, Emil-Fischer-Str. 31, 97074 Würzburg, Germany

45 Department of Physics and Astronomy, University of Turku, Finland, FI-20014 University of Turku, Finland

46 Department of Physics, TU Dortmund University, Otto-Hahn-Str. 4, 44227 Dortmund, Germany

47 INFN Sezione di Roma La Sapienza, P.le Aldo Moro, 2, 00185 Rome, Italy

48 ILANCE, CNRS – University of Tokyo International Research Laboratory, University of Tokyo, 5-1-5 Kashiwa-no-Ha Kashiwa City, Chiba 277-8582, Japan

49 Physics Program, Graduate School of Advanced Science and Engineering, Hiroshima University, 1-3-1 Kagamiyama, Higashi-Hiroshima City, Hiroshima 739-8526, Japan

50 INFN Sezione di Roma Tor Vergata, Via della Ricerca Scientifica 1, 00133 Rome, Italy

51 University of Split, FESB, R. Boškovića 32, 21000 Split, Croatia

52 Department of Physics, Yamagata University, 1-4-12 Kojirakawamachi, Yamagata-shi 990-8560, Japan

53 Institut für Theoretische Physik, Lehrstuhl IV: Plasma-Astroteilchenphysik, Ruhr-Universität Bochum, Universitätsstraße 150, 44801 Bochum, Germany

54 Sendai College, National Institute of Technology, 4-16-1 Ayashi-Chuo, Aoba-ku, Sendai city, Miyagi 989-3128, Japan

55 Josip Juraj Strossmayer University of Osijek, Department of Physics, Trg Ljudevita Gaja 6, 31000 Osijek, Croatia

56 INFN Dipartimento di Scienze Fisiche e Chimiche – Università degli Studi dell’Aquila and Gran Sasso Science Institute, Via Vetoio 1, Viale Crispi 7, 67100 L’Aquila, Italy

57 Chiba University, 1-33, Yayoicho, Inage-ku, Chiba-shi, Chiba 263-8522, Japan

58 Kitashirakawa Oiwakecho, Sakyo Ward, Kyoto 606-8502, Japan

59 FZU – Institute of Physics of the Czech Academy of Sciences, Na Slovance 1999/2, 182 21 Praha 8, Czech Republic

60 Laboratory for High Energy Physics, École Polytechnique Fédérale, CH-1015 Lausanne, Switzerland

61 Astronomical Institute of the Czech Academy of Sciences, Bocni II 1401, 14100 Prague, Czech Republic

62 Faculty of Science, Ibaraki University, 2 Chome-1-1 Bunkyo, Mito, Ibaraki 310-0056, Japan

63 Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies, LPNHE, 4 place Jussieu, 75005 Paris, France

64 Graduate School of Science and Engineering, Saitama University, 255 Simo-Ohkubo, Sakura-ku, Saitama city, Saitama 338-8570, Japan

65 Institut of Particle and Nuclear Studies, KEK (High Energy Accelerator Research Organization), 1-1 Oho, Tsukuba 305-0801, Japan

66 INFN Sezione di Trieste and Università degli Studi di Trieste, Via Valerio 2 I, 34127 Trieste, Italy

67 Escuela Politécnica Superior de Jaén, Universidad de Jaén, Campus Las Lagunillas s/n, Edif. A3, 23071 Jaén, Spain

68 Saha Institute of Nuclear Physics, Sector 1, AF Block, Bidhan Nagar, Bidhannagar, Kolkata, West Bengal 700064, India

69 Institute for Nuclear Research and Nuclear Energy, Bulgarian Academy of Sciences, 72 boul. Tsarigradsko chaussee, 1784 Sofia, Bulgaria

70 Department of Physics and Astronomy, Clemson University, Kinard Lab of Physics, Clemson, SC 29634, USA

71 Institut de Fisica d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra, (Barcelona), Spain

72 Grupo de Electronica, Universidad Complutense de Madrid, Av. Complutense s/n, 28040 Madrid, Spain

73 Institute of Space Sciences (ICE, CSIC), and Institut d'Estudis Espacials de Catalunya (IEEC), and Institució Catalana de Recerca I Estudis Avançats (ICREA), Campus UAB, Carrer de Can Magrans, s/n, 08193 Bellatera, Spain

74 Department of Physics, Konan University, 8-9-1 Okamoto, Higashinada-ku Kobe 658-8501, Japan

75 School of Allied Health Sciences, Kitasato University, Sagamihara, Kanagawa 228-8555, Japan

76 RIKEN, Institute of Physical and Chemical Research, 2-1 Hirosawa, Wako, Saitama 351-0198, Japan

77 Charles University, Institute of Particle and Nuclear Physics, V Holešovičkách 2, 180 00 Prague 8, Czech Republic

78 Division of Physics and Astronomy, Graduate School of Science, Kyoto University, Sakyo-ku, Kyoto 606-8502, Japan

79 Institute for Space-Earth Environmental Research, Nagoya University, Chikusa-ku, Nagoya 464-8601, Japan

80 Kobayashi-Maskawa Institute (KMI) for the Origin of Particles and the Universe, Nagoya University, Chikusa-ku, Nagoya 464-8602, Japan

81 Graduate School of Technology, Industrial and Social Sciences, Tokushima University, 2-1 Minamijosanjima, Tokushima 770-8506, Japan

82 INFN Sezione di Pisa, Edificio C – Polo Fibonacci, Largo Bruno Pontecorvo 3, 56127 Pisa, Italy

83 Gifu University, Faculty of Engineering, 1-1 Yanagido, Gifu 501-1193, Japan

84 Department of Physical Sciences, Aoyama Gakuin University, Fuchinobe, Sagamihara, Kanagawa 252-5258, Japan

85 Macroarea di Scienze MMFFNN, Università di Roma Tor Vergata, Via della Ricerca Scientifica 1, 00133 Rome, Italy