1 Introduction↩︎

Extragalactic sources of TeV \(\gamma\)-rays may be surrounded by lower energy extended emission originating from electrons and positrons deposited in the intergalactic medium between the source and the Earth by interactions of the highest energy \(\gamma\)-rays with low-energy photons of the extragalactic background light [1], [2]. The extended emission appears as an excess in the radial profiles of the source signal [3], [4] that is superimposed on the parent point source signal distributed according to the telescope Point Spread Function (PSF). Searches for such extended emission have so far mostly resulted in the upper limits on the extended emission flux. These upper limits have been used to derive a lower bound on the strength of Intergalactic Magnetic field (IGMF) [5]–[7]. A notable exception is the result of [8] that has reported a detection of the extended emission around the second-brightest extragalactic TeV source on the sky, Mrk 501, a blazar type Active Galactic Nucleus (AGN). This extended emission is found in the data of Fermi Large Area Telescope (LAT). According to [8], such a detection corresponds to a specific combination of IGMF parameters: strength \(B\sim 10^{-15}\) G and correlation length \(\lambda_B\sim 10\) kpc.

Sensitivity of searches for the extended emission with both space-based Fermi/LAT and ground-based Imaging Atmospheric Cherenkov Telescopes [5], [8]–[12] is limited by the precision of knowledge of the telescope PSF. Imperfections in the early modeling of the Fermi/LAT PSF have previously led to false positive detection(s) of extended emission around blazars [13]. At the initial stages of operation of Fermi/LAT telescope, the PSF modeling was based on Monte-Carlo simulations of \(\gamma\)-ray events that did not properly reproduce the observed event distributions. [14] has reported a phenomenological study of Fermi/LAT PSF based on real data, that has revealed limitations of the Monte-Carlo modeling of the PSF. This study has concluded that the detection of extended emission claimed in [13] was spurious. The Fermi/LAT collaboration has performed an independent analysis of the PSF based on real telescope data [15] and reached a similar conclusion. Most recently, the Fermi/LAT collaboration has re-derived the properties of the PSF from improved Monte-Carlo modeling¹ so that the Instrument Response Functions (IRF) released by the Fermi/LAT collaboration currently do not rely on the real data. The Fermi/LAT performance web-page² specifies that a validation study of the PSF parameterisations has been performed and involved a comparison of the PSF models derived from the IRF parameters with the real data of Vela pulsar and of a set of AGN. This web page quotes a 5% level uncertainty on the 68% containment radius of the PSF in the 0.1 GeV to 10 GeV energy range, but no other tests. Useful as it is (see e.g. [15]), the information on the 68% signal containment radius is not sufficient for assessment of the PSF uncertainties in the context of the search of the extended emission around extragalactic sources and discovery of IGMF. It is possible to find two model PSF profiles that have the same 68% containment radii but are largely different in shape (e.g. a Gaussian and an exponential function, or a profile with a power-law tail at large angular distances).

In what follows we compare the Fermi/LAT PSF given by the IRF with observational data with the goal to assess the precision of the PSF parameters knowledge. In section 2 we summarize the approach for description of the PSF in the Fermi/LAT IRF. In section 3 we present our method of measurement of PSF uncertainties using the signals of bright pulsars. We show that the PSF models with parameters given in the IRF files do not properly describe the data. In Section 4 we present an improved PSF that is consistent with pulsar data. We use the improved PSF in re-analysis of Mrk 501 data in section 5. Finally in section 6 we discuss implications of our results.

2 Fermi/LAT PSF parameterisation↩︎

Fermi-LAT IRFs use an analytical model of the PSF to describe the radial distributions of Monte-Carlo simulated events around sources. The PSF is modeled as a sum of two King functions, \[K(x,\sigma,\gamma)=\frac{1}{2\pi\sigma^2}\left(1-\frac{1}{\gamma}\right)\left[1+\frac{1}{2\gamma}\frac{x^2}{\sigma^2}\right]^{-\gamma} \label{eq:king}\tag{1}\] with different parameters \(\sigma_c,\gamma_c\) for the "core" and \(\sigma_t,\gamma_t\) for the "tail" components of the PSF: \[PSF(x)=f_{c}K(x,\sigma_{c},\gamma_{c})+(1-f_{c})K(x,\sigma_{t},\gamma_{t}) \label{eq:psf}\tag{2}\] The factor \(f_c\) describes the relative importance of the core and tail components in the overall PSF.

The parameters \(\sigma_c,\gamma_c,f_c,\sigma_t,\gamma_t\) depend on the energy \(E\) and off-axis angle \(\theta\) of gamma-ray events detected by LAT. These parameters are stored as lookup tables (on a grid of \(E\) and \(\cos\theta\) values) in the IRF files for different event selections. The parameter \(x=\psi/S_P(E)\) in (2 ) is a re-scaled angular distance \(\psi\) from the source with the scaling factor \(S_P(E)\) that depends on energy (and on the event selection) \[S_P(E)=\sqrt{c_0^2\left(\frac{E}{100MeV}\right)^{-2\beta}+c_1^2}\] The re-scaling parameters \(c_0, c_1,\beta\) are also stored in the IRF files. The PSF parameters for different \(\gamma\)-ray event types are stored in different IRF files. Overall, the PSF parameterisations are presented as lookup tables for a set of 23 logarithmically spaced energy bins between 10 MeV and 3 TeV, 8 linearly spaced bins of cosine of the off-axis angle between 0.2 and 1. The \(\gamma\)-ray events are divided into groups based on either type of interaction in the LAT detector (FRONT / BACK) or on the quality of reconstruction of the arrival direction (PSF0-PSF3). For the PSF0-PSF3 grouping, there is the total of 805 parameters, if one considers also the 69 parameters of the function \(S_P(E)\) in each energy bin. For the FRONT/BACK division, there are 437 parameters in total.

To compute the PSF for a specific source in a given sky direction \((RA, Dec)\) at a given energy \(E\), one needs to sum the King functions (1 ) over different off-axis angles, with weights proportional to the expected number of events from the source with different off-axis angles. This number of events as a function of the off-axis angle can be estimated based on the history of telescope pointings (the "spacecraft file" in Fermi/LAT data terminology). Such weighting would produce correct PSF prediction for a source that is not variable in time. Otherwise, a bursting source for which most of the event statistics are accumulated during a short outburst may have a PSF that is very different from the steady-state source PSF, because most of the source signal would be composed of events incident at a specific off-axis angle at which the source was visible during the outburst. This suggests a different approach for determination of the PSF for a bursting source (this approach would also work for a steady-state source). The off-axis angle weighting in the PSF calculation can be done based on the real statistics of events from the source, for example, the events within certain limiting angular distance from the source position.

In general, the uncertainties of the PSF parameters include statistical and systematic uncertainties. The statistical uncertainties are determined by the statistics of the Monte-Carlo event sample and they can in principle be reduced via increase of the Monte-Carlo event sample size. To the contrary, the systematic uncertainties cannot be assessed and cannot be reduced using the Monte-Carlo data alone. There may be certain unaccounted instrumental effects that are not considered in Monte-Carlo simulations, the telescope hardware may have characteristics that evolve with time in a way that is not completely controlled etc. Hence, the systematic uncertainty of the PSF has to be addressed via comparison with real data. Such a comparison should ideally be based on the data on genuine point sources detected by Fermi/LAT.

3 PSF study with pulsars↩︎

The point source signals are naturally provided by pulsars, fast rotating neutron stars. In general, pulsar signals are superimposed onto extended emission from its pulsar wind nebula and on diffuse Galactic emission. There may also be weak sources in the direct vicinity of the pulsar and it may be challenging to deconvolve emission from these weak sources from the pulsar emission. In this sense, pulsar total emission does not represent a "clean" point source. However, the extended pulsar wind nebula flux, the weak nearby point sources and the diffuse Galactic emission signals do not vary periodically and are identical during the "on-pulse" and "off-pulse" phases of the pulsed signal. The only difference in the signal between the two phases is the presence of the flux coming from the direct vicinity of the neutron star, within the distance of the "light cylinder" which is clearly not resolvable by Fermi/LAT for any of the known pulsars. Subtracting the "off-pulse" flux from the "on-pulse" flux automatically removes the pulsar wind nebula, the weak nearby sources and the diffuse Galactic emission. It is well possible that the "off-pulse" signal also contains weaker emission from the pulsar magnetosphere and subtraction of the "off-pulse" flux can also remove part of the point source flux. This reduction of the point source flux reduces the statistics of the "pure" point source signal, but does not affect its morphology. Thus the off-pulse subtracted signal from pulsars provides a genuine point source flux that can be used for study of the systematic uncertainty of the PSF.

We use pulsars to study the PSF models provided in Fermi/LAT IRFs. We do this for the specific event class P8R3_SOURCEVETO_V3. This event selection is characterized by the lowest residual charged cosmic ray background and is therefore most suitable for the studies of weak extended sources (as specified in recommendations on data selections by Fermi Science Support Center, https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/ lat_data_selection.html). The characteristics of the Monte-Carlo derived PSF for this event class are given in the psf_P8R3_SOURCEVETO_V3_PSF.fits (for the PSF0-PSF3 event grouping) and psf_P8R3_SOURCEVETO_V3_FB.fits (for the FRONT-BACK event grouping) files of Fermi/LAT IRF set for the Pass 8, Version 3 data (P8V3).

Our analysis considers stacked signal of the three brightest pulsars on the sky, Vela, Geminga and Crab. We select events of different event types, PSF0,...,PSF3, FRONT, BACK, by filtering LAT event files using gtselect-gtmktime tools chain (see https://fermi.gsfc.nasa.gov/ssc/data/analysis/) specifying the event types that correspond to different event selections. We retain only events in the energy range above 100 MeV. We limit our data-PSF model comparison to the regions within \(5^\circ\) around the pulsar position. We consider events up to maximum Zenith angle \(100^\circ\). At the gtmktime step, we impose a filter (DATA_QUAL)&&(LAT_CONFIG==1), as recommended by the Fermi Science Support Center (https://fermi.gsfc.nasa.gov/ssc). We also limit the time span of the data for each pulsar according to the time limits of validity of available pulsar ephemerides [16].

To select the "on-pulse" and "off-pulse" signals, we use pulsar ephemerides provided with the 3-rd pulsar catalog of Fermi/LAT [16]. We use the pint program [17] to assign the pulse phase for each event. The ephemerides files are taken from the data repository associated with the pulsar catalog paper (https://fermi.gsfc.nasa.gov/ssc/data/access/lat/ 3rd_PSR_catalog/). The current version of pint is not able to handle the ephemerides file for the Crab pulsar from this repository and we use the tcb2tdb script to convert the time format of the ephemerides file from "TCB" to "TDB". Plotting the event data as a function of the pulse phase we obtain the pulse profiles as shown in Fig. 1 for the energy bin 1-1.7 GeV. Based on these pulse profiles, we define the "on-pulse" and "off-pulse" phases. For Vela, we use the on-pulse phase interval \(0.05<\phi<0.675\), for Geminga \(-0.05<\phi<0.675\) and for Crab \(0.025<\phi<0.575\). The phase interval complementary to the on-pulse is the off-pulse interval used for estimation of the background on top of which the point source signal is observed.

As an example, Fig. 2 illustrates the difference in the event types PSF0,...,PSF3 for different quality of the angular reconstruction. The two panels show events binned in equal angular bins \(\psi\) and equal solid angle bins \(\psi^2\) respectively. The top panel gives a better insight on the behavior of the inner PSF while the bottom panel provides a view of the "surface brightness" profile of the PSF. One can see that the events belonging to the PSF3 type have the best angular reconstruction quality and are distributed closest to the source. The events from the type PSF0 have the worst PSF quality. The angular width of the PSF for all the event types increases with decreasing energy and below 100 MeV it becomes larger or comparable to the angular cut of \(5^\circ\) considered in our analysis. We limit the analysis to \(E>100\) MeV for this reason (only bins starting from energy bin 4 of the IRF file are considered). The pulsed emission signal statistics becomes weak in the energy range above 10 GeV and our analysis is effectively limited to the energy range below this energy.

Fig. 3 and 4 show comparisons of the model PSF radial profiles from Fermi/LAT IRFs with the pulsar data in the energy ranges 1-1.7 GeV and 5.6-10 GeV for different event types. From the top panels of the figures one can see that the analytical model described in section 2 provides an overall satisfactory description of the data, in the first approximation. The bottom panels show residuals of the data fits with the PSF models derived from the IRF file, with positive residuals shown in black and negative residuals in red.

These panels show discrepancies of models with the data, that may reach way beyond 10% of the PSF over a wide range of angles, as for example in the last panel of Fig. 4. Even though the fraction of the point source signal contained in the PSF tails is small, the data – model mismatch cannot be ignored in the studies of extended emission, because the mismatch mimics an excess over the point source signals that can be produced by the extended emission.

Deviations of the shapes of the model PSF from the real data can be spotted in the analysis of residuals of the model fits to the pulsar data. Fitting the PSF models to the data for each energy bin, each off-axis angle and each event type, we find that most of the time, the data are not consistent with the model of the PSF given by the IRF and better model fits can be found by fitting the PSF model parameters to the pulsar data.

Fig. 5 shows the difference in the \(\chi^2\) between the data fits with the PSF models with parameters extracted from the IRF file and with the five PSF parameters \(f_c,\sigma_c,\gamma_c,\sigma_t,\gamma_t\) left free in the fitting. With improvements of \(\chi^2\) by as much as \(\Delta\chi^2\ge 100\) in more than half of the energy and off-axis angle bins, the PSF parameters given in the IRF file lie outside the \(90\%\) confidence range of the \(f_c,\sigma_c,\gamma_c,\sigma_t,\gamma_t\) measurements from the pulsar data. The white contours in Fig. 5 show these 90% contours that correspond to \(\Delta\chi^2=10.7\) given that the fit adjusts six parameters: the five parameters of the PSF plus the overall normalization of the model [18]. The mismatch of the PSF model parameters and the data in individual off-axis angle bins is illustrated in Fig. 6.

4 Improved PSF model based on pulsar data↩︎

Limitations of the PSF models derived from Monte-Carlo studies, described in the previous section, motivate an alternative data-driven approach for modeling of the PSF in the energy range 0.1-10 GeV where high statistics signal from pulsars is available.

One possibility is to use directly the pulsed emission signal, rather than model PSF, as a template for the point source signals. We bin the real pulsar events into the same energy and off-axis angle bins as the IRF file describing the PSF. Data points in Figs. 6, 7 show examples of the binned pulsar signals for specific energy, off-axis angle and PSF type bins. One can see that the statistics of the pulsar signal at 100 MeV – 1 GeV energy range is so high that the error bars are barely noticeable in the profile within \(1^\circ-2^\circ\) distance from the source, seen if the data are split into narrow off-axis angle bins. Replacing the analytical models of the PSF with the real-data PSF has an advantage that it completely removes the systematic error of the PSF modeling, including possible error introduced through the choice of the analytical template of the PSF: King function vs. Gaussian, one King function vs. two- or three- King function etc. models.

Using the binned pulsed emission data as PSF templates, it is possible to construct the expected PSF for a source in a specific sky direction, based on the observed distribution of events within certain distance from the source as a function of the off-axis angle, \(W(\cos\theta)\). The PSF in a given energy bin is found by summing the lookup table \(T\) over the off-axis angle \[\label{eq:psf95convolution} PSF(E,\psi)=\sum_{\cos\theta}W(\cos\theta)T_{pulsar}(E,\cos\theta,\psi)\tag{3}\]

An alternative possibility is to re-use and improve the analytical model of the PSF described in section 2 by fitting the parameters \(f_c,\sigma_c,\gamma_c,\sigma_t,\gamma_t\) to the pulsed emission data. An example of such fitting is shown in Fig. 6, 7 for a specific energy and off-axis angle bins by solid lines surrounded by gray shadings that show the 68% uncertainty ranges of the model fits. To calculate these error bands, we have run Monte-Carlo simulations randomly displacing the PSF parameters from the best fit and noticing the PSF models for which the \(\chi^2\) of the model fit does not exceed the best-fit model \(\chi^2\) by more than \(\Delta\chi^2=7.04\) (which corresponds to the \(68\%\) confidence range for the six fitted parameters [18]). These improved PSF models are compared with the initial PSF models from the IRF. Different types of the model-data discrepancies can be spotted for the PSF models from the IRF: a mismatch in the core of the PSF is visible in Fig. 6, a mismatch in the PSF tail is seen in Fig. 7. There are no such discrepancies in the re-fitted models that match the data by definition.

As mentioned in the Introduction, PSF models can differ in shape (as shown in Fig. 7 and Fig. 6) but can have the same 68% containment radii, meaning that the containment radius is possibly not the only parameter useful for the model validation. This is illustrated by Fig. 8 from which one can see that the 68% containment radii for both PSF agree with each other up to few percent.

Table ¿tbl:tab:psf95lookup? shows an example comparison of the original and revised PSF parameter tables for a specific energy bin, \(1-1.7\) GeV, for the FRONT/BACK event grouping. The table also gives the \(\chi^2\) of the pulsar data fits with the model PSF from the IRF file and best-fit model.

An advantage of the analytical PSF model compared to the stacked data PSF model is that it provides a smaller error at large angular distances from the source. This is clearly seen from Figs. 6, 7: the gray uncertainty bands of the analytical models are narrower than the scatter of the stacked signal data points at large angular distances from the source. A disadvantage of the analytical PSF models is the presence of a residual systematic error due to the fact that the analytical model does not necessarily provide the perfect description of the PSF shape. For the rest of our study, we consider the analytical PSF model only.

5 Application of revised PSF to the search of extended emission around Mrk 501↩︎

Mrk 501 is the second-brightest blazar on the TeV sky [19]. Detectability of extended secondary emission around this source has been studied in [20], which concluded that magnetic field as strong as \(10^{-11}\) G is measurable through the search of the extended secondary emission in the energy range up to \(E>100\) GeV. The range of weaker magnetic field strengths can be explored via the search of the secondary signal at lower energies accessible to Fermi/LAT, the idea that has been exploited in [8]. It is worth noticing that Mrk 501 is not the best source for the search of the extended emission in the energy range of Fermi/LAT, because it has a very strong point source signal in 1-100 GeV range. Sources with harder spectra, like the extreme blazars 1ES 0229+200 or 1ES 0502+675 provide better sensitivity because of lower ratio between the expected point source and extended emission fluxes [6].

The secondary flux at the energy \(E_\gamma\sim 10\) GeV is produced following interactions of \(\gamma\)-rays with energies \[E_{\gamma_0}\simeq 3\left[\frac{E_\gamma}{10 GeV}\right]^{1/2} TeV\] in the intergalactic medium [3]. The mean free path of such \(\gamma\)-rays is \[D_{\gamma_0}\sim 2\times 10^2\left[\frac{E_{\gamma_0}}{3 TeV}\right]^{-0.5} Mpc\] in the range \(1 TeV<E_{\gamma_0}<10 TeV\) [21], [22]. This is larger than the distance to the source at the redshift \(z\simeq 0.034\). The secondary emission detectable in Fermi/LAT energy range is accumulated all over the Line-of-Sight (LoS) toward Mrk 501.

If the IGMF has parameters \(B\sim 10^{-15}\) G and \(\lambda_B\sim 10\) kpc, as derived in [8], the angular size of the extended source detectable with Fermi/LAT in 1-10 GeV range is expected to scale with energy as [3] \[\Psi_{ext}=1.2^\circ\left[\frac{B}{10^{-15} G}\right]\left[\frac{\lambda_B}{10 kpc}\right]^{1/2}\left[\frac{E_\gamma}{10 GeV}\right]^{-3/4} \label{eq:theta95ext}\tag{4}\] which is comparable to the extent of the Fermi/LAT PSF that has the 95% containment radius \(\Psi_{PSF}\simeq 1^\circ\) at 10 GeV (see https://www.slac.stanford.edu/exp/glast/groups/canda/ lat_Performance.htm) for the event selection used in [8].

We start our analysis of Mrk 501 by reproducing the results of [8] with the original PSF parameterisations with parameters extracted from the IRF files, using fermipy (v1.4.0)³ [23] and fermitools (v2.2.0)⁴ packages. We subsequently replace the PSF from the original IRF files with the updated PSF based on our pulsar analysis and redo the same analysis using the revised PSF. For the energies outside the pulsar data range, we consider the original PSF parameters.

We select Fermi-LAT data from 17-years of operation between March 2008 and January 2026 in the energy range 0.1-300 GeV. We define a \(10^{\circ}\times 10^{\circ}\) ROI centered on Mrk 501 with spatial bin size of \(0.02^{\circ}\) and 4 bins per energy decade. We follow the recommended data selection of the Fermi-LAT website for extragalactic diffuse analysis of Pass 8 Data (P8R3)⁵ and use the P8R3_SOURCEVETO_V3 IRF (\(\texttt{evclass=2048}\)) with FRONT+BACK photons (\(\texttt{evtype=3}\)) filtered with the corresponding time selection filtering (DATA_QUAL\(>\)​0)&&(LAT_CONFIG==1). We run the gtselect-gtbin-gtmktime-gtexpcube2-gtscrmaps tool chain to initialize the analysis setup. Finally, we use the best-fit synchrotron-self Compton emission model of [8] for the spectrum of Mrk 501 defined as \[\label{eq:Mrk50195int95spec} \frac{dN}{dE} = N_0 \left(\frac{E}{E_0}\right)^{-\Gamma}\exp{\left(-\frac{E}{E_c} \right)^{\beta}},\tag{5}\] with \(N_0=1.04 \times 10^{-5}~\mathrm{TeV^{-1}s^{-1}cm^{-2}}\), \(\Gamma=1.72\), \(E_c=0.33~\mathrm{TeV}\), and \(\beta=0.37\).

We perform a likelihood ratio test where we compute the test statistic \[\label{eq:TS} TS=-2(\ln L_0-\ln L_1)~,\tag{6}\] where \(\ln L_0\) and \(\ln L_1\) are the log-likelihoods to be maximized in the following two hypotheses.

In the null hypothesis \(H_0\), we consider the sky model that contains all sources from the 4FGL-DR4 Fermi catalog [24] within a \(15^{\circ}\times 15^{\circ}\) Region-Of-Interest (ROI) centered on Mrk 501. We free the normalization and spectral model parameters of all sources included in a radius of \(5^{\circ}\) around Mrk 501 or with a \(TS\geq25\). We additionally free the normalization of the isotropic and galactic diffuse models and fix all spectral model parameters of Mrk 501 except its normalization. We fit the ROI model by maximizing the log-likelihood.

In the alternative hypothesis \(H_1\), we use the same initial configuration as in \(H_0\) and add a source representing the extended emission around Mrk 501 simulated with CRbeam code [25] for a source with the intrinsic spectrum given by Eq. (5 ). The primary gamma-ray source is assumed to emit photons within a jet of \(\theta_{\text{jet}}=3^{\circ}\) opening angle. The jet direction is aligned with the line of sight. Secondary photons with a time-delay below \(10^7\) years are binned in spatial and energy bins to yield a source template that is used in the likelihood analysis. We allow the normalization of the extended source \(\epsilon>0\) to be free in the fit. \(\epsilon=1\) corresponds to the secondary flux predicted by the model for specific magnetic field parameters which we take to be the best-fit values of the analysis of [8], \(B=1.5\times 10^{-15}\) G and \(\lambda_B=10\) kpc. The best-fit \(\epsilon\) parameter is rescaled according to the best-fit normalization of Mrk 501.

The analysis based on the default PSF from the IRF files gives \(\ln L_0=-861849.18\) and \(\ln L_1 = -861838.72\) translating to \(TS=20.92\) with \(\epsilon=0.89\pm 0.12\), which corresponds to a \(\simeq4.57\sigma\) significance of detection of the extended emission. This is somewhat lower than the significance reported in [8]. The difference can be due to a difference in the modeling of the extended emission (we use CRbeam code, while the analysis of [8] was based on ELMAG code [26]) and due to the difference in the event selection: SOURCEVETO photons vs. SOURCE photons.

Once we replace the standard PSF with the PSF derived from the pulsar analysis, we obtain \(\ln L_0=-861895.21\) and \(\ln L_1 = -861893.17\) corresponding to the best-fit \(\epsilon=0.40\pm 0.12\) with the test statistic \(TS=4.08\), so that the extended emission is not significantly detected anymore. To the contrary, the TS value of the point source at the location of Mrk 501 increases from \(TS=49375.54\) to \(TS=49579.00\), which hints for a better match between the revised PSF to the point source signal shape.

The maximal likelihood of the fit with the ROI model based on revised PSF is slightly below that of the model based on the IRF PSF. To check if the model based on the revised PSF is consistent with the data, we have performed the "goodness-of-fit" test, by running a set of Monte-Carlo simulations to produce mock Fermi-LAT data realizing the sky model, using simulate_roi function of fermipy package. Fitting the mock data with the reference model, we find a distribution of the log-likelihood values expected in the case when the model is consistent with the data. This distribution is shown in Fig. 9. We see that the log-likelihood found for the real data is well within the expected range.

To infer an upper limit on the extended secondary source flux (at 95% confidence level, CL), we look at a deviation of \(\Delta TS=-2.71\) from the \(H_0\) hypothesis (\(TS=0\)) by scanning the \(\epsilon\) parameter space and re-fitting the ROI model at each step. The results are shown on the left panel of Fig. 10 where the maximal halo normalization for the revised analytical PSF model is \(\epsilon_{\text{upper}}=0.93\), which is lower than the expected flux normalization \(\epsilon=1\). Hence, the extended emission for the best-fit values \(B=1.5\times 10^{-15}\), \(\lambda_B=10\) kpc of the IGMF is ruled out at the \(95\%\) CL. Additionally, we present the results for the default analytical PSF model on the right panel of Fig. 10 where \(\epsilon_{\text{upper}}=2.00\), indicating no disagreement with the presence of the extended emission model from CRbeam code at the 95% CL. Fig. 11 shows the 95% CL upper limit on the extended emission flux derived from Fig. 10 with the default and revised PSF models supporting the same conclusions.

6 Discussion and conclusions↩︎

Knowledge of the telescope PSF is crucial for the search of extended \(\gamma\)-ray emission around extragalactic sources that may be induced by the effects of propagation of \(\gamma\)-rays through the intergalactic medium. Use of imprecise PSF model may lead to appearance of false positive detection of extended emission that would be favored by a fit of observational data to "patch" the mismatch between the data and the PSF model in certain angular distance range. Past examples of false positive detections of extended emission around blazars in Fermi/LAT data [13] have previously led to revision of Fermi/LAT PSF models driven by the model-data comparisons [14], [15].

In this paper we performed the model-data comparison for the PSF models provided by the IRFs of Pass 8, revision 3 event selection of Fermi/LAT data. We have found discrepancies between the models and data on pulsed emission from bright pulsars in the energy range 0.1-10 GeV. The model parameters quoted in the IRF files happen to be inconsistent with the data in different energy and off-axis angle bins. We have revised the PSF modeling for a specific event class, P8R3_SOURCEVETO, that is most suitable for the search of the extended emission around extragalactic sources.

We have demonstrated the importance of the precise PSF modeling for the search of extended emission around gamma-ray sources by re-assessing the result on detection of extended emission around Mrk 501 [8]. Our analysis shows that replacing the PSF models derived from the analytical parameterisations provided by the IRFs with the data-driven PSF "washes out" an evidence for the presence of the extended emission around this source, previously claimed in [8].

The data-driven approach for the PSF determination has a limitation related to the necessity of high-statistics genuine point source signal provided by pulsars. Such signal is only available in the energy range below \(\sim 10\) GeV. Above this energy, the pulsed emission fluxes of pulsars are suppressed and the statistics of the signal drops dramatically. This means that our method for PSF characterization would not work at the highest energies and the Monte-Carlo based PSF determination would still be needed for the extended emission searches up to the highest energies accessible to Fermi/LAT.

References↩︎