Geminga

General information

• Name of the source: Geminga
• Brief description of the source: Radio-quiet middle-aged gamma-ray pulsar
  • Object type : Pulsar
  • Other relevant information: Third pulsar detected by IACTs
  • RA: 06 33 54.15 (hh mm ss), Dec: +17 46 12.9 (dd mm ss)
  • RA, Dec in deg (ICRS): 98.48 17.77
  • Galactic coordinates: l=195.13 b=+04.27

People involved

• Álvaro Mas-Aguilar (alvmas@ucm.es)
• Marcos López-Moya (marcos@gae.ucm.es)
• Giulia Brunelli (giulia.brunelli@inaf.it)
• Rubén López-Coto (rlopezcoto@iaa.es)
• Giovanni Ceribella (ceribell@mpp.mpg.de)
• Paul K. H. Yeung (pkh91yg@icrr.u-tokyo.ac.jp)

Data-taking information

Geminga has been observed with standard wobble observations (offset = 0.4 deg), in dark and good weather conditions:

• TNG dust < 0.5 ug/m3
• No clouds both from satellite images and all sky camera
• Rate_LL around 5.5 kHz
• Lidar trans@9km >= 0.85

Runs per each day

• 2022

 20-12-2022: [11533 11534 11535 11536 11537 11538 11539 11540 11541]
 22-12-2022: [11580 11581 11582 11583 11584 11585]

• 2023

 14-01-2023: [11645 11646 11647 11648]
 15-01-2023: [11664 11665 11666]
 17-01-2023: [11693 11694 11695 11696 11697 11698]
 18-01-2023: [11715]
 20-01-2023: [11732 11733 11734 11735 11736 11737]
 21-01-2023: [11761 11762]
 24-01-2023: [11818 11819 11820 11821 11822 11823 11824 11825]
 25-01-2023: [11850 11851 11852 11853 11854 11855]
 26-01-2023: [11884 11885 11886 11887]
 28-01-2023: [11950 11951 11952]
 11-02-2023: [11959 11960]
 13-02-2023: [11972 11973 11974 11975 11976]
 14-02-2023: [11983 11984 11985 11986 11987 11988 11989]
 18-02-2023: [11998]
 20-02-2023: [12018 12019 12020 12021 12022 12023 12024]
 22-02-2023: [12041 12042 12043 12044]
 24-02-2023: [12070 12071 12072 12073 12074]
 10-03-2023: [12139 12141]
 13-03-2023: [12152 12153 12154 12155 12156 12157 12158 12159 12160]
 15-03-2023: [12178 12179 12180 12181]
 16-03-2023: [12197 12198 12199 12200 12201]
 18-03-2023: [12248 12249 12250 12251 12252 12253]
 19-03-2023: [12283 12284 12285 12286 12287 12288]
 21-03-2023: [12340 12341 12342 12343 12344]
 22-03-2023: [12373 12374 12375 12376 12377]
 23-03-2023: [12385 12386]
 24-03-2023: [12415 12416 12417 12418 12419 12420 12421]
 07-12-2023: [15867 15868 15869 15870 15871 15872]
 08-12-2023: [15899 15900 15901 15902 15903 15904 15905 15906 15907]
 09-12-2023: [15931 15932 15933 15934 15935 15936 15937 15938 15939]
 10-12-2023: [15960 15961 15962 15963 15964 15965 15966 15967 15968]
 12-12-2023: [16029 16030 16031 16032 16033 16034 16035 16036 16037]
 13-12-2023: [16067 16068]
 14-12-2023: [16100 16101 16102 16103 16104 16105 16106 16107 16108 16109]
 15-12-2023: [16138 16139 16140 16141 16142 16143 16144 16145 16146]
 17-12-2023: [16210 16211 16212 16213 16214 16215]
 18-12-2023: [16232 16233 16234 16235]
 19-12-2023: [16252 16253 16254]

• 2024

 09-01-2024: [16282 16283 16284]
 10-01-2024: [16310 16311]
 14-01-2024: [16370 16371 16372 16373 16374 16375 16376]
 01-02-2024: [16459 16460 16461 16462 16463 16464 16465 16466 16467 16468]
 04-02-2024: [16549 16550 16551 16552 16553 16554 16555]
 05-02-2024: [16578 16579 16580 16581]
 06-02-2024: [16584 16585 16586 16588 16589]
 11-02-2024: [16700 16701 16702 16703 16704 16705 16706 16707 16708 16709]

Total observational time without any selection: 79.25h

The quality selection was based on the differential intensity spectra of the detected showers, dR/dI (events/s/p.e.), see reference notebook here. A zenith cut Zd<50deg was also applied.

Final list of selected runs:

[11533, 11534, 11535, 11536, 11537, 11538, 11539, 11580, 11581, 11582, 11583, 11584, 11585, 11645, 11646, 11647, 11648, 11664, 11665, 11666,
 11693, 11694, 11695, 11696, 11697, 11698, 11715, 11732, 11733, 11734, 11735, 11736, 11737, 11761, 11762, 11818, 11819, 11820, 11821, 11822,
 11823, 11824, 11825, 11850, 11851, 11852, 11853, 11854, 11855, 11884, 11885, 11886, 11887, 11950, 11959, 11960, 11972, 11973, 11974, 11976,
 11983, 11984, 11985, 11986, 11987, 11988, 11989, 12018, 12041, 12042, 12043, 12044, 12139, 12141, 12152, 12153, 12154, 12155, 12156, 12157,
 12158, 12159, 12197, 12198, 12199, 12200, 12201, 12248, 12249, 12250, 12251, 12252, 12283, 12284, 12285, 12286, 12287, 12288, 12340, 12341,
 12342, 12343, 15867, 15868, 15869, 15870, 15871, 15872, 15899, 15900, 15901, 15902, 15903, 15904, 15905, 15906, 15907, 15931, 15932, 15933,
 15934, 15935, 15936, 15937, 15938, 15961, 15962, 15963, 15964, 15965, 15966, 15967, 15968, 16029, 16030, 16031, 16032, 16033, 16034, 16035,
 16036, 16037, 16067, 16068, 16100, 16101, 16102, 16103, 16104, 16105, 16106, 16107, 16108, 16109, 16138, 16139, 16140, 16141, 16142, 16143,
 16144, 16145, 16146, 16210, 16211, 16212, 16213, 16214, 16215, 16232, 16233, 16234, 16235, 16282, 16283, 16284, 16459, 16460, 16461, 16462,
 16463, 16465, 16466, 16467, 16468, 16549, 16550, 16551, 16552, 16553, 16554, 16579, 16580, 16584, 16586, 16588, 16589, 16700, 16701, 16702,
 16703, 16704, 16705, 16706, 16707, 16708]

Total good-quality observational time: 60.1h

Data analysis

Low-level analysis

Software version used: lstchain v0.10.7

Approach: source-dependent

DL1 to DL2

The analysis started from the DL1a files produced by LSTOSA:

/fefs/aswg/data/real/DL1/20*/v0.10/tailcut84/dl1*.h5

A dedicated source-dependent production of Monte Carlo data was required (see here). The MCs were produces at dec=22.76deg and with NSB tuned on Geminga. The corresponding RF models used for the DL2 production were stored in:

/fefs/aswg/data/models/AllSky/20240126_v0.10.4_src6_dec2276_tuned/dec_2276/

Path to the DL2 files:

/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/src_dependent_v0.10/DL2/ 

DL2 to IRF

The corresponding MC Testing Dataset was used for the IRF production:

/fefs/aswg/data/mc/DL2/AllSky/20240126_v0.10.4_src6_dec2276_tuned/TestingDataset/dec_2276/

The point-like IRFs were computed using a linear interpolation, with the following cuts

• gh-efficiency=70%
• alpha-containment=70%
• intensity > 50 p.e.

specified in the configuration file

/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/src_dependent_v0.10/irf_tool_config.json

Path to the IRFs:

/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/src_dependent_v0.10/IRFs/

DL2 to DL3

The same configuration file was used for the DL3 production.

Path to the DL3 files:

/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/src_dependent_v0.10/DL3/

DL3 to Pulsar DL3

To produce the final DL3 including the phase information, the software PulsarTimingAnalysis (available here) has been used. The software uses the package PINT (version v0.9.7) to compute the phase information, using a user-provided ephemeris file, and adds it to the DL3 files.

The ephemeris used in the analysis were provided by Giovanni Ceribella and are available here. The ones used for the analysis were updated on 31/03/2024 and can be found here:

/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/src_dependent_v0.10/ephemeris_updated_31032024.par

Path to the Pulsar DL3 files:

/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/src_dependent_v0.10/DL3_pulsar/

High-level analysis

Software version used: gammapy v1.1

Phaseogram

To produce the phaseogram of Geminga, the software PulsarTimingAnalysis, based on gammapy, has been used.

The phase regions selected for the two peaks, the bridge and the background are:

P1: [0.056, 0.161]
P2: [0.550, 0.642]
Bridge: [0.161, 0.550]
OFF: [0.700, 0.950]

Binning:

• 60 bins in phase for all phaseograms
• 5 bins per decade for the phaseograms per energy range

P2 morphology

To study the morphology of P2, two symmetric profiles, the Gaussian and the Lorentzian, were tested. They are described by

fG(φ)=AG exp[ - (φ - μ)2/2σ2] + constG

for the Gaussian and

fL(φ)=AL [ 1 + (φ - x0 / 2γ)2]-1 + constL

for the Lorentzian.

We fitted with scipy.optimize.curve_fit() both the LST-1 and the Fermi-LAT phaseograms in two cases:

• without applying any energy cut (i.e. full-band fit)
• in the reconstructed energy range [15, 65] GeV, where most of the signal of P2 is detected, dividing into two logarithmically-spaced bins, [15,31] GeV and [31, 65] GeV.

For each fit, we derived the mean pulse position (μ and x₀) and the Full-Width-Half-Maximum (FWHM), computed as

FWHMG=2σ sqrt(2 ln2)

for the Gaussian and

FWHMG=2γ

for the Lorentzian, and we computed the Akaike Information Criterion (AIC = 2k − 2 log(L)) to assess the preference for either one of the two models.

Spectral fit of P2

The fit has been done only for P2, stacking the LST-1 data.

Energy binning:

e_reco = MapAxis.from_energy_bounds(0.01, 10, nbin=40,unit="TeV", name="energy")
e_true = MapAxis.from_energy_bounds(0.003, 50, nbin=100, per_decade=False, unit="TeV", name="energy_true")

We employed a mask in the spectral fit based on a subset of the e_reco axis, whose lower and upper edges are 19.95 GeV and 94.5 GeV. An extension of the same subset (16.7-94.5 GeV) has been used to estimate the Flux points.

We employed a Power Law spectral model

dN/dE=N0(E/E0)-Γ

for the main fit, but also a Log-Parabola model

dN/dE=N0(E/E0)-a-blog(E/E0)

to test for the presence of curvature in the energy range.

Spectral fit of P1

We attempted at obtaining some flux points or upper limits on the P1 spectrum. To do so, we fitted with a Power Law in the same energy range as P2 and 3 energy bins. We fixed the spectral index to Γ=4.5 and the reference energy to 25 GeV, otherwise the fit does not converge.

This test was done just to compare the estimated upper limits to the Fermi-LAT flux points, but it is not certain whether P1 is a fluctuation or not.

Joint fit with Fermi-LAT data

A new analysis of 14 years of Fermi-LAT data was performed to obtain the phaseograms and the flux points. Then, the flux points were employed to perform a joint fit with the stacked LST-1 spectral points (starting at the DL3 level) using gammapy.

The Fermi flux points are stored in the files:

/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/Geminga_P2_sed.csv

for P2 and

/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/Geminga_P1_sed.csv

for P1.

Systematics estimation

For the evaluation of the systematics, four effects have been considered:

• Change of the MC efficiency
• Introduction of an energy scaling factor in the IRFs
• Change of the Zd cut
• Change of the intensity cut

Change of the MC efficiency

Both the alpha containment and the gh efficiency have been changed, considering the values of:

• gh-eff = 40% - 70% - 90%
• alpha-cont = 70% - 90%

For each combination of gh-eff and alpha-cont, the reduction from the DL2 and IRF stage has been performed. All the data products (DL3, IRFs and Pulsar DL3), as well as the configuration files, are available inside the corresponding directories:

/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/src_dependent_v0.10/ghX_alphaY/

where X is the gh-eff and Y is the alpha-cont.

Introduction of an energy scaling factor in the IRFs

An energy scaling of ±15% has been added when creating the IRFs and the DL3 files. To do so, a test branch of lstchain has been used, for which the configuration file had an additional parameter scale_true_energy. Setting "scale_true_energy": 0.85 corresponds to the -15% scaling, while "scale_true_energy": 1.15 corresponds to the +15% scaling.

The reduction has been run from the IRF production stage and all the data products are available in the folders:

/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/src_dependent_v0.10/systematics_irfs/en_scale_minus15/
/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/src_dependent_v0.10/systematics_irfs/en_scale_15/

Change of the Zd cut

This effect has been studied at the high-level stage of the analysis, changing the cut in zenith when building the list of the selected runs to produce the spectrum on gammapy.

The following zenith cuts have been tested: 50°, 40°, 30°, 25°, 20°, and 15°.

Change of the intensity cut

The intensity cut of the data has been changed, testing values between 20 p.e. and 70 p.e., and performing the reduction from the IRF level.

All the data products, and respective configuration files, are available here:

/fefs/aswg/workspace/giulia.brunelli/data/real/Geminga/src_dependent_v0.10/systematics_intensity/intX/

where X is the value of the cut.

Results

Phaseogram

Phaseogram obtained for the complete 2022-2024 sample at Zd<50°, without any energy cut (Figure 1 of the paper).

Phaseograms for different energy bins (Zd<50°)

P2 spectrum

Spectrum of Geminga P2 in the range [20, 95] GeV.

P2 morphology

Best-fit results, with the associated statistical error, for the P2 mean position and width (FWHM) for the symmetric Gaussian and Lorentzian profiles. We report the results for both LST-1 and Fermi-LAT, in the two energy bins and the full band. The Aikake Information Criterion (AIC) for every fit is shown. 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 (Table 1 of the paper).

P2 joint fit with Fermi-LAT

Joint LST-1 and Fermi-LAT data samples of P2, along with the best-fit results of both the power law with an exponential cutoff (PLEC, dotted line) and the power law with sub-exponential cutoff (PLSEC, dashed line). The power law fit of the LST-1 only points (orange squares) is shown together with its statistical error bar (solid line and shaded area). The MAGIC Collaboration+ (2020) points are reported as green triangles as a comparison (Figure 3 of the paper).

Systematic estimation

The values of the index and of the normalisation obtained for each systematic effect are reported in the following tables:

• change of the MC efficiency

• addition of energy scaling

• change of the zenith distance

• change of the intensity cut

Comparison with Harding+ (2021) model

SED points of P2 obtained from the LST-1 sample, along with its best-fit power law model and its associated statistical uncertainty band, overlayed to the Synchro-Curvature component of the Harding et al. (2021) model for Geminga. Fermi-LAT and MAGIC points from the analysis presented in this work and in MAGIC Collaboration et al. (2020), respectively, are also reported for comparison. The solid black line corresponds to the model with values of the low and high accelerating electric field of R^low_acc = eE_low/mc² = 0.04 and R^high_acc = eE^high/mc²=0.15, while the dashed line is the model with R^low_acc = R^high_acc = 0.15. The plot has been adapted from Harding et al. (2021), where more details can be found. (Figure 4 of the paper)

Attempt of P1 spectrum

Presentations at internal meetings

• 30/01/2023, first report on the detection (LST-reco call, PDF)
• 19/05/2023, update on the analysis of the 2022-2023 sample (LST-gal call, PDF)
• 22/05/2023, update on analysis of 2022-2023 sample (LST-reco call, PDF)
• 19/06/2023, update on the analysis of the 2022-2023 sample (LST GM Munich, PDF)
• 09/11/2023, update on the analysis of the 2022-2023 sample (LST GM Fall 2023, PDF)
• 20/11/2023, cross-check with source-independent analysis of the 2022-2023 sample (LST-reco call, PDF)
• 29/01/2024, update on the cross-check with source-independent analysis of the 2022-2023 sample (LST-reco call, PDF)
• 18/03/2024, update on the source-dependent analysis of the total sample (LST-reco call, PDF)
• 22/03/2024, update on the source-dependent analysis of the total sample (LST-gal call, PDF)
• 21/05/2024, update on the analysis of the total sample (LST GM Prague, PDF)
• 10/06/2024, update on the systematics estimation (LST-reco call, PDF)
• 17/06/2024, update on the systematics estimation (LST-reco call, PDF)
• 27/09/2024, presentation of the final results (LST-gal call, PDF)
• 30/09/2024, presentation of the final results (LST-reco call, PDF)
• 13/11/2024, presentation of the final results and paper updates (LST GM Fall 2024, PDF)
• 21/11/2024, presentation of the final results and paper updates (CTAO Consortium Galactic WG call, PDF)

Presentations at external conferences/meetings

Presentation