© The Authors 2024

1 Introduction

The Fermi Large Area Telescope (LAT) is by far the most sensitive γ-ray telescope in the giga-electronvolt band ever flown (Atwood et al. 2009). Since its launch in 2008, it has continuously monitored the γ-ray sky and detected thousands of new sources. The latest source catalog (4FGL-DR4; Ballet et al. 2023) is based on 14 yr of data and supersedes three previous iterations of the fourth Fermi-LAT source catalog (Abdollahi et al. 2020, 2022; Ballet et al. 2020). It contains 7195 γ-ray emitters from 50 MeV to 1 TeV, distributed over the whole sky. The two most prominent contributing source classes are blazar-type active galactic nuclei (AGN) and rotation-powered pulsars (PSRs), of which around 3900 and 320 can be found in the catalog, respectively. However, a large number of sources, around 2600, do not yet have a convincing association at lower energies, providing considerable discovery potential in the Fermi-LAT catalogs. In particular, the distribution of unassociated Fermi sources in galactic latitude suggests a likely Galactic origin for a significant fraction of the 2600 unidentified sources. While a significant fraction of these low-latitude sources may originate from mismodeled diffuse emission (Abdollahi et al. 2022), it is a promising effort to search for new pulsars in individual unassociated Fermi sources, as numerous previous studies have shown (e.g., Au et al. 2023; Bruzewski et al. 2023; Hare et al. 2022; Braglia et al. 2020).

The sample of known γ-ray pulsars consists of approximately equal numbers of young, energetic (i.e., non-recycled) pulsars and millisecond (i.e., recycled) pulsars, all of whose emission is powered by rotation rather than accretion (Smith et al. 2023). A considerable fraction of the young energetic pulsars is found to be radio-quiet (Abdo et al. 2013; Smith et al. 2023), so that their γ-radiation is often necessary for their detection. On the other hand, a large fraction of MSPs are found in binary systems with an optical companion, due to their formation requiring spin-up by accretion in a low-mass X-ray binary (Alpar et al. 1982; Radhakrishnan & Srinivasan 1982). An especially interesting class of MSP binaries is given by black-widow and redback systems (see Roberts 2013; Hui & Li 2019), which are characterized by a strong interaction between the pulsar and its companion and short orbital periods. In these “spider” systems, strong particle winds from the pulsar continually ablate their companion star (Fruchter et al. 1988; Strader et al. 2019), often leading to enhanced X-ray emission from intra-binary shocks (Gentile et al. 2014; Koljonen & Linares 2023).

It is well established that rotation-powered pulsars typically emit X-rays at a level of $L_{\text{X}}/\dot{E}\sim {10}^{-5}-{10}^{-3}$ of their spin-down power (Becker & Trümper 1997; Shibata et al. 2016), with frequent contributions from thermal emission from the neutron star surface and nonthermal magnetospheric emission. Typical giga-electronvolt γ-ray efficiencies are significantly larger, at $L_{\text{X}}/\dot{E}\sim {10}^{-2}-1$ (Abdo et al. 2013). For Fermi-detected pulsars in particular, the giga-electronvolt γ-ray and X-ray fluxes are correlated with large scatter, likely due to a variety of emission mechanisms, viewing geometries, and beaming fractions (Marelli et al. 2011, 2015). Typical flux ratios (or equivalently luminosity ratios) lie in the range F_γ/F_X ~ 100–10 000, with young radio-quiet pulsars being relatively X-ray fainter than radio-loud ones and MSPs (Abdo et al. 2013). In contrast to radio, X-ray, or γ-ray energies, pulsars exhibit comparatively weak emission at optical or infrared wavelengths (Abdo et al. 2013; Mignani 2011).

The biggest obstacle in identifying counterparts to Fermi-LAT sources is their inherently large positional error. This makes a direct search for counterparts, for instance, in the optical, challenging, as the expected number of positional chance coincidences is enormous. Therefore, a frequently employed approach in the search for possible associations for Fermi sources is the usage of X-ray data, as a large fraction of blazars and rotation-powered pulsars are expected to show detectable X-ray emission. The positions of X-ray sources are usually accurate on an arc-second level, facilitating associations with counterparts at radio or optical wavelengths. An excellent example of this is the systematic follow-up project of unassociated Fermi sources with Swift-XRT by Stroh & Falcone (2013)¹. Within this program, X-ray observations of Fermi error ellipses have been routinely carried out to search for potential counterparts, leading to many newly detected pulsars. While snapshot surveys such as that by Stroh & Falcone (2013) provide immense discovery potential, they are limited by the field of view of their telescope, which prevents the reliable estimation of the probability of chance coincidences.

The eROSITA telescope aboard the Spektrum-Roentgen-Gamma (SRG) observatory (Predehl et al. 2021; Sunyaev et al. 2021) provides the first all-sky survey in X-rays since the ROSAT mission (Truemper 1982) around 30 yr prior. Within its bandpass of 0.2–10 keV, the eROSITA all-sky survey (eRASS) is providing the deepest coverage of the X-ray sky to date, with an expected final (i.e., after the eight planned surveys) all-sky point-source sensitivity in the 0.2–2.3 keV band around 25 times deeper than the ROSAT all-sky survey (Predehl et al. 2021; Merloni et al. 2012). At the time of writing, four out of eight planned all-sky surveys have been completed, with the cumulative dataset (named eRASS:4) reaching a typical point source sensitivity of around 1–2 × 10⁻¹⁴ erg s⁻¹ cm⁻² in the main band. This paper employs the eRASS:4 data of the western Galactic hemisphere, that is, in the range 180° < l < 360°. This dataset provides the deepest contiguous picture of half the X-ray sky, including the regions around 1400 unassociated Fermi sources, with almost three million detected X-ray sources.

The main goal of this work is to identify the X-ray counterparts of (candidate) rotation-powered pulsars in eRASS:4 data. On the one hand, this is achieved via a probabilistic cross-match of eROSITA X-ray sources to unassociated Fermi sources, which are likely to be high-energy pulsars. This is achieved by combining a purely positional cross-match with a prior classification of unassociated 4FGL sources based on their γ-ray properties, and constraints on the expected X-ray flux of a potential match. With this effort, we aim to provide useful lists of promising targets for radio, optical, X-ray, or γ-ray follow-up by quantifying the likelihood of an association with a given X-ray source using a Bayesian approach (Budavári & Szalay 2008; Salvato et al. 2018). On the other hand, we present the set of detected X-ray counterparts of known rotation-powered pulsars in eRASS:4, emphasizing those sources whose X-ray emission was observed for the first time to extend the sample of X-ray-detected pulsars.

Our paper is organized as follows: we give an overview of the data products we used in Sect. 2, and we describe the individual components of our method for constructing probabilistic cross-match catalogs in Sect. 3. In Sect. 4, we present the sample of X-ray counterparts to previously identified pulsars to verify our capabilities of detecting X-ray emission from pulsars. In Sect. 5, we present our list of candidate pulsar matches to unassociated γ-ray sources and discuss their implications and follow-up prospects. Finally, we summarize our work in Sect. 6.

2 Input data from Fermi-LAT and eROSITA

The target of our analysis is the set of all unidentified sources in the Fermi-LAT 14-yr source catalog (4FGL-DR4, in the following 4FGL), which contains 7195 giga-electronvolt γ-ray sources (Ballet et al. 2023). From this catalog, we selected all unidentified sources, including those with low-probability associations, meaning all sources labeled in the catalog with CLASS1 == “”, “unk”. In order to characterize the properties of the typical γ-emitting source classes, we also defined two samples of associated sources, a pulsar sample (CLASS1 == “psr”, “msp”), and a sample of blazar-type AGN (CLASS1 == “bll”, “fsrq”, “bcu”). This results in a total of 2577 unassociated sources, of which 1371 lie within the western Galactic hemisphere. Among the complete set of associated sources, 320 are of pulsar- and 3934 of blazar-type according to our definition, of which 163 and 1922 are located in the western Galactic hemisphere, respectively.

The primary eROSITA data product we use is the latest internal eRASS:4 single-band source catalog². This catalog contains almost three million X-ray sources with a 0.2–2.3 keV detection likelihood (Brunner et al. 2022) above 5. The catalog columns important for this work are the sky positions (RA_CORR, DEC_CORR), and positional error (RADEC_ERR) of the sources, as well as the estimated X-ray flux and its error (ML_FLUX_1, ML_FLUX_ERR_1). Furthermore, for the purpose of visual inspection of X-ray sources, we used the available eRASS:4 imaging data in the regions of all Fermi sources in the latest processing version c020. We manually rebinned the available images from the archive to create simple count images with a pixel size of 12″ for each region of interest. The energy bands we used for the image creation were the three standard eROSITA bands 0.2–0.6, 0.6–2.3, and 2.3–5.0 keV.

3 Cross-matching strategy

Typically, the extent of 4FGL error ellipses is very large (the median 95% semi-major axis being 5.8′ for unassociated sources), and the density of X-ray sources in the all-sky survey is quite high (at an average source density per sky area of 140 deg⁻² in eRASS:4). Hence, one expects to observe many chance coincidences between 4FGL and eROSITA sources in a naive positional cross-match, with a median of three X-ray sources expected per 4FGL error ellipse. We, therefore, followed a more complex scheme, ranking potential associations between the two catalogs in a quantitative, probabilistic manner. Apart from a positional cross-match (Sect. 3.2), our approach included a prior classification of 4FGL sources based on their γ-ray properties (Sect. 3.1), obtaining constraints on the expected X-ray flux of a pulsar/blazar-type match (Sect. 3.3), and the visual inspection of the resulting matches and possible multiwave-length counterparts (Sect. 3.5). Finally, we tested our method using the known associations of 4FGL sources (Sect. 3.6).

3.1 Classification of γ-ray sources

It has been shown numerous times (Salvetti 2014; Saz Parkinson et al. 2016; Kaur et al. 2019; Kerby et al. 2O2la,b; Germani et al. 2021; Bhat & Malyshev 2022; Coronado-Blázquez 2022), that the properties of γ-ray sources in the Fermi-LAT source catalogs can be used to predict their source type, and in particular to differentiate between pulsar-like and blazar-like sources. Typically, the primary differences between the two classes were found to be a stronger spectral downward curvature for pulsars, and stronger variability for blazars. Here, we performed a machine-learning-based classification of unassociated 4FGL-DR4 sources, with the main target of estimating a prior probability $P_{\gamma }^{\text{PSR}}$ for a particular source to be a pulsar, rather than an AGN. While this neglects the presence of further source classes among the sample of unassociated 4FGL sources, we believe this approach is sufficient to separate pulsar candidates from the most prominent expected background population, as only 8% of associated 4FGL sources are assigned to such further classes.

The type of classifier we used for this task was a random forest classifier, as was frequently done before (Salvetti 2014; Saz Parkinson et al. 2016; Bhat & Malyshev 2022). In order to tackle the strongly unbalanced nature of the training set (with a factor of ten more blazars than pulsars), we decided to use a BalancedRandomForest classifier (Chen et al. 2004), as implemented in the imbalanced-learn package (Lemaître et al. 2017). Briefly, this tool combats the class imbalance by undersampling the majority class, rather than artificially creating synthetic data points to oversample the minority class (e.g., SMOTE, Chawla et al. 2002).

The process we followed to tune our classifier is explained in detail in Appendix A. Importantly, we attempted to avoid introducing biases against faint sources, by excluding absolute fluxes from our list of candidate features and by constructing flux-independent versions of certain features (e.g., curvature significance, flux error). In order to determine the optimal set of features for our classifier, we performed recursive feature elimination (Luo et al. 2020), with the target of maximizing its average precision (see Appendix A), meaning the ability of the classifier to construct a pure and complete sample of pulsar candidates from a test sample. The optimal set of features, along with their associated importance, is displayed in Table 1, and the corresponding precision-recall curve, determined via repeated five-fold cross-validation, is depicted in Fig. 1. Our classifier achieved an average precision of 0.920, and a ROC AUC³ score of 0.989, which is at a similar level as similar classifiers based on earlier source catalogs. For comparison, the work of Saz Parkinson et al. (2016) reports a ROC AUC of 0.993.

The result of our classifier is a binary probability $P_{\gamma }^{\text{PSR}}$ for each unassociated source, quantifying the likelihood of it being a pulsar rather than an AGN based on its γ-ray properties. The distribution of this quantity is compared to that obtained for pulsar and AGN test samples in Fig. 2. Overall, our classifier yields a relatively large fraction of sources with $P_{\gamma }^{\text{PSR}}>0.5$, namely around 45%, compared to other works (e.g., Bhat & Malyshev 2022; Coronado-Blázquez 2022). This is likely caused by our approach of using a balanced classifier, which implicitly assumes that both classes are a priori equally likely, as well by our attempt to not bias ourselves against faint sources (see Appendix A). Figure 3 shows the distribution of pulsar or blazar candidates identified by our classifier across the Galaxy. As expected, a large fraction of sources in the Galactic plane are classified as likely pulsars, whereas the majority of sources off the plane are found to be likely AGN. Nonetheless, there is a noteworthy presence of high-latitude pulsar candidates, which are particularly promising for identification due to less crowded fields at any wavelength.

In order to classify whether individual pulsar candidates are likely to be young or MSPs, we analogously tuned and trained a second classifier on the corresponding pulsar subsamples in the 4FGL-DR4 catalog (141 young pulsars vs. 179 MSPs). This classifier (see Table A.2 for the used features) returns a probability $P_{\gamma }^{\text{YNG}}$ for a given source to be young (rather than recycled), conditional on its pulsar nature. This allows us to identify, for instance, likely MSP candidates in our sample, which could constitute promising sources for optical follow-up. Our algorithm achieved a ROC AUC score of 0.964 on the test sample during 5-fold cross-validation, preforming similarly to previous attempts (Saz Parkinson et al. 2016).

All results of our γ-ray based prior classification efforts are provided in electronic form at CDS. In the supplied table, for each unassociated 4FGL source, we give the raw parameters used by our classifiers, and the resulting binary classification probabilities $P_{\gamma }^{\text{PSR}}$ and $P_{\gamma }^{\text{YNG}}$, as defined in this section.

Fig. 1 Precision-recall curve for the binary classifier separating pulsar-and blazar-candidate γ-ray sources. This graph compares the recall of the classifier, i.e., the fraction of pulsars in the test sample that are correctly identified, to its precision, i.e., the fraction of true pulsars among all objects selected as pulsars, dependent on the decision threshold.

Fig. 2 Distribution of binary classification probabilities $P_{\gamma }^{\text{PSR}}$ for the targeted unassociated 4FGL sources. This is compared to the probabilities evaluated on PSR- and AGN-type “out-of-box” test samples during cross-validation, meaning targets the classifier did not see during training. The fractions of misclassified AGN and pulsars in the test sample are 4.2% and 4.9%, respectively.

3.2 Positional cross-match

For the probabilistic quantification of positional agreement between two sources in the eRASS:4 and 4FGL-DR4 catalogs, we follow an approach equivalent to that outlined by Budavári & Szalay (2008) and implemented in NWAY (Salvato et al. 2018):

The prior probability for the hypothesis H that a randomly selected eROSITA source is the correct counterpart of a given 4FGL source is given by $$P(H)=\frac{c}{v\text{Ω}},$$(1)

where v describes the sky density of sources in the eROSITA catalog, Ω refers to the overlapping sky area covered by the two catalogs, and c is the prior completeness factor, describing the fraction of Fermi sources for which we expect to find an X-ray counterpart (Salvato et al. 2018). Given the observed data D, the likelihood for a match of two sources with a positional separation ψ, and (Gaussian) positional errors σ₁ and σ₂ is given by (see Eqs. (16) and (28) in Budavári & Szalay 2008): $$P(D\mid H)=\left(\frac{\text{Ω}}{4\pi }\right)\times \frac{2}{{\sigma }_1^2+{\sigma }_2^2}\exp \left(-\frac{{\psi }^2}{2\left({\sigma }_1^2+{\sigma }_2^2\right)}\right).$$(2)

Hence, following from Bayes’ theorem, the posterior probability of a given match is proportional to $$P(H\mid D)\propto \frac{c}{2\pi v\left({\sigma }_1^2+{\sigma }_2^2\right)}\exp \left(-\frac{{\psi }^2}{2\left({\sigma }_1^2+{\sigma }_2^2\right)}\right).$$(3)

In practice, given that the Fermi-LAT source positions are specified with error ellipses, we used the full covariance matrix Σ to model the astrometric error accurately. Furthermore, given the considerable variation of the X-ray source density on the sky, we replaced the constant sky density of sources with a location-dependent one v = v(x), which we estimated locally within a two-degree radius of each Fermi source. This yields the following expression for the astrometric component of our posterior: $$\begin{array}{c}P(H\mid D)\propto \frac{c}{2\pi v(x)\sqrt{\left|\Sigma \right|}}\exp \left(-\frac{1}{2}{\psi }^{\top }{\Sigma }^{-1}\psi \right),\text{ where }\\ \Sigma =\left(\begin{array}{cc}{\sigma }_{1,x}^2+{\sigma }_2^2& {\rho }_1{\sigma }_{1,x}{\sigma }_{1,y}\\ {\rho }_1{\sigma }_{1,x}{\sigma }_{1,y}& {\sigma }_{1,y}^2+{\sigma }_2^2\end{array}\right).\end{array}$$(4)

Here, we have decomposed the separation vector ψ into its x and y components. The positional covariance matrix components for 4FGL sources (ρ₁, σ_1,x, σ_1,y) were determined from the given 95% error ellipses (rescaled by a factor 0 409 to obtain a “one-sigma”-ellipse) following Pineau et al. (2017). The quantity σ₂ corresponds to the one-sigma positional error of the eRASS:4 sources, determined from their RADEC_ERR following the scheme of Merloni et al. (2024).

3.3 Including X-ray and γ-ray flux information

The majority of all sources in eRASS:4 are detected with only a few counts, the median value of background-subtracted counts being around 10 in the 0.2–2.3 keV band. Therefore, we are unable to aid our classification with detailed analyses of spectral or timing properties of X-ray sources. However, it is well known that rotation-powered pulsars and blazars exhibit vastly different γ-ray to X-ray flux ratios, with blazars being much X-ray brighter in a relative sense (Kaur et al. 2019; Kerby et al. 2021b). Furthermore, phenomenological differences exist also between different types of pulsars, with radio-quiet young pulsars tending to be X-ray fainter than radio-loud ones and MSPs (Marelli et al. 2011, 2015; Abdo et al. 2013). Therefore, depending on the pulsar type (young or MSP), the expected X-ray flux of its potential counterpart may vary drastically.

Hence, with the goal of predicting possible X-ray fluxes for blazars, young pulsars, and MSPs, we constrained the expected flux ratio F_γ/F_X for each class based on the associated 4FGL sources and their X-ray counterparts in our catalog. We cross-matched all counterpart coordinates of associated sources with the eRASS:4 catalog using a 15″ matching radius. For each match, we recorded the X-ray (pseudo-)flux in the 0.2–2.3 keV band (ML_FLUX_1⁴), and γ-ray flux in the 100 Me V–100 Ge V band. In case of no matching source, we used the internal eROSITA upper limit server (see Tubín-Arenas et al. 2024) to determine the counts and background level within a 30″ aperture⁵ at the source position. Cases in which a nearby overlapping source was present were excluded from this analysis, as no reliable upper limits could be extracted here. Overall, we obtained 21, 31, and 1762 usable X-ray fluxes and 42, 46, and 89 upper limits for young pulsars, MSPs, and γ-ray-bright blazars, respectively. The distribution of X-ray and γ-ray flux information for the associations is shown in Fig. 4, and the distribution of X-ray to γ-ray flux ratios in Fig. 5. In addition, we give a list of all X-ray detected counterparts of 4FGL pulsars in Table B.2.

In order to make predictions for fluxes of potential matches of unassociated 4FGL sources, we constructed a simple model of the flux ratio distribution for each source class. We chose to use a log-normal distribution with arbitrary offset and scale for this purpose log F_γ/F_X ~ 𝒩(µ, σ). To determine the best-fit parameters of the distribution, we employed an unbinned likelihood approach for the X-ray detected sources, where the likelihood of a given set of parameters (µ, σ) is given by $$\log ℒ(\mu ,\sigma )={\displaystyle \underset{i}{{\sum }^{\text{}}}}\log \mathcal{N}\left(\log F_{\gamma ,i}/F_{X,i}\mid \mu ,\sigma \right),$$(5)

where the sum runs over all detected sources i. For undetected sources, meaning sources with X-ray upper limits only, we used the full Poissonian likelihood of the intrinsic X-ray flux (see Tubín-Arenas et al. 2024) together with the observed γ-ray flux, to estimate the likelihood of possible flux ratios for a given source. By combining the likelihoods of X-ray detected and undetected sources, we obtained an estimate of the intrinsic flux ratio distribution of each source class, which is unbiased by the eRASS:4 sensitivity limit. For each class, Fig. 5 displays the distribution of observed flux ratios F_γ/F_X and their lower limits for X-ray nondetections, together with the best-fitting log-normal distribution, whose parameters are given in Table 2.

Our analysis qualitatively confirms several previously established findings. This includes the much brighter appearance of blazars compared to pulsars in X-rays, indicated by the much lower flux ratios. In Fig. 5, this can be observed in the much larger fraction of detected AGN compared to young pulsars and MSPs, and also in the evident difference between the peaks of the distributions of detected sources. Furthermore, young pulsars clearly show a much larger spread in F_γ/F_X than MSPs. This may be caused by a significant contribution of very X-ray faint radio-quiet pulsars to the population of young γ-ray pulsars (Abdo et al. 2013). We note, however, that a significant discrepancy may exist between the cataloged X-ray fluxes, which are based on count rates and a simple energy conversion factor (Merloni et al. 2024), and the intrinsic source fluxes, especially in the case of very soft or hard spectra, or strong foreground absorption. In addition, since the main eROSITA energy band is very soft, there is likely a significant contribution of thermal X-ray emission from the neutron star surface, while, for instance, Marelli et al. (2011) only considered its nonthermal component, which may lead to a relative decrease in F_γ/F_X here. Furthermore, the nonthermal X-ray and γ-ray emission of pulsars is likely to be nonisotropic, as significant beaming is expected (Becker 2009). Since the spatial origin of magnetospheric X-ray and γ-ray emission may be different, different beaming fractions and viewing geometries may affect the two regimes, with variations between young and recycled pulsars being possible (Marelli et al. 2011). All these factors may introduce offsets and scatter into the observed F_γ/F_X relation, compared to what may be the relation of their intrinsic magnetospheric luminosities.

In order to evaluate whether a given candidate match is a likely X-ray counterpart to an unassociated 4FGL source, it is necessary to compare our expected class-dependent flux distribution with the flux distribution of field X-ray sources in the region of interest. To achieve this, for each γ-ray source, we constructed a histogram of X-ray source fluxes within a 2° radius and fit a simple empirical model for the background flux distribution ϕ_B of the type: $${\varphi }_{\text{B}}\left(\log F_{\text{X}}\mid A,\text{Γ},\lambda ,\theta \right)=A{\left(F_{\text{X}}\right)}^{1+\text{Γ}}\text{Φ}\left(\log F_{\text{X}}\mid \lambda ,\theta \right),$$(6)

where Φ denotes a cumulative normal distribution, which models the sensitivity-limited cutoff (as in Ge et al. 2022). Our fit yields constraints on the characteristic slope of the local X-ray flux distribution Γ, on the logarithm of the limiting X-ray flux at which half the sources are detected λ, and on the characteristic width of this cutoff θ. The constant A normalizes the probability distribution to unity. In a few, mostly Galactic, cases, we found that a single power law could not describe the intrinsic distribution of X-ray fluxes well. Hence, for each region, we systematically tested whether adding a second power law or a log-normal component could improve the fit, keeping the model with the lowest Akaike information criterion (Akaike 1974) in each case.

By applying this sensitivity cutoff to the expected class-dependent X-ray flux distribution for a given γ-ray flux F_γ, we can determine the expected probability distribution function ϕ_t of observing an X-ray counterpart to a source of class t at a given flux F_X: $${\varphi }_t\left(\log F_{\text{X}}\mid \log F_{\gamma }\right)=\text{Φ}\left(\log F_{\text{X}}\mid \lambda ,\theta \right)\mathcal{N}\left(\log F_{\text{X}}/F_{\gamma }\mid {\mu }_t,{\sigma }_t\right).$$(7)

Hence, by evaluating the quotient ϕ_t(log F_X | log F_γ)/ϕ_B(log F_X) for any potential match, we obtain a factor modifying its likelihood (see Budavári & Szalay 2008; Salvato et al. 2018), under the assumption of a certain source class. Furthermore, integrating the function ϕ_t over all possible X-ray fluxes yields an estimate of the prior completeness, that is, the probability of a given γ-ray source having a detected X-ray counterpart (see Sect. 3.4). The approach of comparing the expected flux distributions for counterparts and field sources is illustrated for an example 4FGL source in Fig. 6.

This step of constraining the expected flux of the X-ray counterpart to a given Fermi source is vital to avoid overestimating the likelihood of having detected a pulsar-type match. X-ray counterparts to γ-ray pulsars frequently have a small chance of having a flux detectable by eRASS:4, in particular for fainter γ-ray sources (F_γ ≲ 3 × 10⁻¹² erg s⁻¹ cm⁻²). The converse is true for blazars, which are typically expected to stand out as a particularly bright X-ray source compared to the field.

To test whether rudimentary spectral information of our X-ray sources could be used to aid our classification, we also inspected the band-to-band hardness ratios for the energy ranges 0.2–0.5, 0.5–1.0, 1.0–2.0, and 2.0–5.0 keV for our different source types. We did find some global differences between the populations, with pulsar counterparts exhibiting larger, and γ-ray blazars smaller, scatter than field X-ray sources unrelated to the 4FGL catalog. However, due to the weak nature of these trends in our sample, we were unable to construct a quantitatively reliable predictor estimating the likelihood of a particular source being a pulsar- or AGN-type counterpart to a 4FGL source.

Fig. 3 Distribution of binary classification probabilities $P_{\gamma }^{\text{PSR}}$, indicated by the color scale of the markers, of unassociated 4FGL sources across the Galaxy. The gray shading marks the eastern Galactic hemisphere, for which no eROSITA X-ray data are available.

Fig. 4 Distribution of X-ray (0.2–2.3 keV) and γ-ray (100 MeV – 100 GeV) fluxes of young pulsars (blue circles) and MSPs (green). One-sigma X-ray upper limits are indicated with triangles in cyan (young pulsars), lime green (MSPs), and magenta (AGN), respectively. In order to provide a clean plot, the distribution of detected AGN is illustrated by the red filled contours, marking its one, two, and three-sigma outlines.

Fig. 5 Distribution of the γ-to-X-ray flux ratio of 4FGL young pulsars (top), MSPs (middle), and AGN (bottom). In both panels, the distribution of the flux ratio of detected sources (dotted histogram) and corresponding one-sigma lower limits (dashed histogram) are indicated, together with a best-fitting log-normal distribution (solid line).

Fig. 6 Example for the computation of the expected X-ray flux for a specific γ-ray source, 4FGL J0159.8–2234. The top panel indicates the normalized X-ray flux distributions of field sources (black), and for the assumption of a young pulsar, MSP, or AGN nature. The bottom panel depicts the corresponding ratio of probabilities. The markers indicate the fluxes of the detected X-ray sources around the Fermi source. For reference, the γ-ray flux of the source is log Fγ = −12.00, and the limiting X-ray flux is λ = −13.95.

3.4 Computation of posterior match probabilities

To “piece together” the information obtained in the previous sections, we calculated a combined posterior probability for a given X-ray source to be the correct match to an unidentified Fermi source and for that source to be of a certain class. Considering all possibly matching X-ray sources and object types, the normalized posterior probability for an individual match is given by: $$\begin{array}{l}{P}_i^t=\frac{{\rho }_i^t}{\overline{c}+\underset{i}{{{\displaystyle \sum }}^{\text{}}}\underset{t}{{{\displaystyle \sum }}^{\text{}}}{\rho }_i^t},\text{ with }\hfill \\ {\rho }_i^t=P_{\gamma }^t\frac{{\varphi }_t\left(\log F_{\text{X}}\mid \log F_{\gamma }\right)}{{\varphi }_{\text{B}}\left(\log F_{\text{X}}\right)}\frac{\exp \left(-\frac{1}{2}{\psi }^{\top }{\Sigma }^{-1}\psi \right)}{2\pi v(x)\sqrt{\left|\Sigma \right|}}\hfill \\ \overline{c}=1-{\displaystyle \underset{t}{{\sum }^{\text{}}}}{P}_{\gamma }^t{{\displaystyle \int }}^{\text{}}{\varphi }_{\text{t}}\left(\log F_{\text{X}}\mid \log F_{\gamma }\right)\text{d\hspace{0.17em}}\log F_{\text{X}},\hfill \end{array}$$(8)

where the index i runs over all possibly associated X-ray sources, and t = {YNG, MSP, AGN} indicates the considered object types. Here, ${\rho }_i^t$ can be understood as an unnormalized match probability, taking into account the γ-ray properties (via the source-type prior $P_{\gamma }^t$), X-ray flux (via the expected flux distribution for an association ϕ_t and the background ϕ_B), and the astrometry of the involved sources. The quantity $\overline{c}$ measures the “prior incompleteness” (Salvato et al. 2018) for a given γ-ray source, that is, the prior probability for its true counterpart to be missing, given its expected X-ray flux and the sensitivity limit of the catalog. For a given γ-ray source, the probability of having any match of a given type t is given by $P_{\text{any }}^t=\underset{i}{{{{\displaystyle \sum }}^{\text{}}}^{\text{}}}{P}_i^t$, and of having any match of any type by $P_{\text{any }}=\underset{i}{{{\displaystyle \sum }}^{\text{}}}\underset{t}{{{\displaystyle \sum }}^{\text{}}}{P}_i^t$.

After having computed the “pulsar-match” probability $P_i^{\text{PSR}}=P_i^{\text{YNG}}+P_i^{\text{MSP}}$ for each possible X-ray association of a 4FGL source, we compiled our final candidate list for counterparts to possible γ-ray pulsars by applying the following cuts:

• Location within “three sigma” from the γ-ray position, that is within the error ellipse containing 99.7% of the probability mass (around 1.4 times the radius of the 95% ellipse in the 4FGL catalog).

• Minimum pulsar-match probability $P_i^{\text{PSR}}\ge 0.02$.

This left us with around 1100 candidate matches, 200 of which have a probability $P_i^{\text{PSR}}\ge 0.1$.

3.5 Visual inspection and multiwavelength counterparts

As a final step of our procedure, we performed a visual inspection of all potentially matching X-ray sources, fulfilling the criteria specified in the previous section. To achieve this, we created diagnostic plots for each matched 4FGL source, overlaying the X-ray and γ-ray source position(s) on local eRASS:4 X-ray images in the energy bands 0.2–0.6, 0.6–2.3, and 2.3–5.0 keV. The target of this effort was to recognize spurious X-ray sources, which were falsely detected, for instance, due to an enhanced level of diffuse emission in Galactic extended sources, or in the wings of the point spread function of very bright point sources (see also Merloni et al. 2024). While such spurious emission is often classified as extended in the detection pipeline, we refrained from generally excluding extended sources from our cross-matching effort, as young pulsars can potentially exhibit X-ray-bright pulsar wind nebulae, which may show a significant angular size.

In addition to examining X-ray images, we also found it helpful to check whether our potentially matching X-ray sources could be classified securely as non-pulsar via their multiwave-length counterparts. The weak optical and infrared emission of pulsars contrasts the properties of stars and AGN, the two most common X-ray-emitting source classes in the all-sky survey (Salvato et al. 2022), which tend to be characterized by bright optical (Schneider et al. 2022) and mid-infrared (Salvato et al. 2018) emission, respectively. This can be exploited to classify the potential counterpart of an X-ray source (e.g., Tranin et al. 2022). Therefore, we overplotted the positions of our X-ray sources with the distribution of sources in the Gaia DR3 (Gaia Collaboration 2016, 2023) and CatWISE2020 catalogs (Marocco et al. 2021), to check for evident non-pulsar counterparts. In this way, purely coronal X-ray emitters are identifiable as comparatively bright optical sources, with significant proper motion or parallax in the Gaia catalog, whereas AGN tend to appear as mid-infrared sources with typical WISE color W1 – W2 > 0 (Salvato et al. 2018), without significant proper motion or parallax in their possible Gaia counterpart. For completeness, we also queried the SIMBAD database (Wenger et al. 2000) to check for existing identifications of astronomical sources in the vicinity. Furthermore, we queried a subset (SUMSS, TGSS, FIRST, NVSS, VLSSR, WENSS, GLEAM, LOTSS, WISH; Mauch et al. 2003; Intema et al. 2017; Helfand et al. 2015; Condon et al. 1998; Lane et al. 2014; Rengelink et al. 1997; Hurley-Walker et al. 2017; Shimwell et al. 2022; De Breuck et al. 2002) of the radio source catalogs used by Bruzewski et al. (2023) to check whether our sources exhibit radio emission at frequencies ≲ 1 GHz, which may be expected for both blazars and pulsars.

We note that there are two fundamental caveats to this approach: Since, due to their formation, most MSPs are expected to be located in binary systems, it is likely that a large fraction of X-ray-emitting rotation-powered pulsars in fact has an optical counterpart detectable by Gaia. However, since coronal X-ray emitters typically emit only a fraction of their luminosity at X-ray energies (Pizzolato et al. 2003; Wright et al. 2011), we expect such cases of optically detected binary companions to be much fainter optically than a coronal X-ray emitter would be. Hence, even for convincing positional matches, we only excluded sources with optical counterparts with a Gaia magnitude of m_G ≲ 15. For comparison, typically, redback companions have optical magnitudes on the order of m_G ~ 20 (Clark et al. 2023). For potential AGN counterparts, we did not apply a limiting magnitude, but required a convincing extragalactic nature (i.e., no significant proper motion or parallax in Gaia), and a clear mid-infrared counterpart. Second, for most sources, the X-ray astrometric errors and counterpart source densities are so high that a secure statement on the presence of a multiwavelength counterpart is impossible to make. Hence, the above criteria were applied rather conservatively, and sources were only excluded if the multiwavelength association was deemed very likely from a positional standpoint. In practice, this implies positional agreement within 1–2 sigma and sufficiently small error bars for only a single source to be within ~3 sigma of the X-ray source.

The diagnostic plots displaying X-ray images of all sources and multiwavelength catalogs around their positions are available on Zenodo⁶. For illustration, Fig. 7 displays an example plot, showing the region of the source 4FGL J1544.2–2554. In that particular example, an optical counterpart appears to exist for the matched X-ray source. While it is much too faint to explain the observed X-ray source with coronal emission, it could be a viable binary counterpart to a possible X-ray bright MSP.

3.6 Verifying our approach with previously associated 4FGL sources

To test the reliability of our method and verify its probabilistic predictions, we also ran our algorithm on samples of previously associated sources in the 4FGL catalog. In two separate cross-matching runs, we applied our method to all previously identified or associated pulsar- and blazar-type γ-ray sources with CLASS1 == ’psr’, ’msp’ and CLASS1 == ’bll’, ’fsrq’, ’bcu’, respectively. Importantly, this approach provides us with a ground truth to test our method against, as the positions of the true multiwavelength counterparts of 4FGL sources are given in the catalog. To emulate the behavior of our machine-learning classifier on unassociated sources, we predicted the γ-ray-based likelihood for a source to be pulsar $P_{\gamma }^{\text{PSR}}$, using our out-of-bag estimates created via 5-fold cross-validation (Sect. 3.1). This means predictions for a particular source were made by a classifier that had not previously seen that source. Furthermore, we naturally used the full γ-ray-only error ellipses to model match probabilities rather than the known counterpart positions. By binning the resulting cross-matches by match probability $P_i^{\text{PSR}}$ or $P_i^{\text{AGN}}$, and evaluating the number of correctly and incorrectly identified counterparts in each bin, we are thus able to quantitatively test the performance and reliability of our algorithm. We considered a given X-ray source matched correctly if its position was less than 15″ away from the counterpart given in 4FGL-DR4. Since pulsars frequently exhibit extended, not necessarily symmetric, pulsar wind nebulae, we used a more generous radius of 60″ for extended X-ray sources.

The result of this exercise is displayed in Fig. 8. In the cases of both pulsars and blazars, the sample of sources with high match probabilities is found to be very pure, as desired. Furthermore, the quantitative reliability of our output probabilities is illustrated by the close correspondence between predicted match probabilities $P_i^{\text{PSR}}/P_i^{\text{AGN}}$, and the observed fraction of correctly identified counterparts ƒ_det, despite a larger scatter for the pulsar class due to the much smaller test sample.

4 Newly discovered X-ray counterparts to known rotation-powered pulsars

In addition to the identification work done on the 4FGL catalog, we correlated the ATNF pulsar database (Manchester et al. 2005) in version 1.71 with the eRASS:4 source catalog, as quite a significant number of pulsars found by Fermi were subsequently detected in the radio band. An earlier correlation by Prinz & Becker (2015) based on archival XMM-Newton and Chandra data identified the X-ray counterpart of 20 previously undetected rotation-powered pulsars. eROSITA with its unlimited field of view in survey mode, allows for searching for X-ray emission from all rotation-powered pulsars down to a homogeneous sensitivity limit. A statistical analysis of source positions in eRASS:l showed that the astrometric accuracy varies between ~4–15 arc-sec for sources detected with a likelihood below 10 (Merloni et al. 2024). We thus restricted the search radius to 15 arcsec, and, using the combined eRASS:4 dataset, we found that 97 pulsars from the ATNF pulsar database in the western Galactic hemisphere match by position with an eROSITA source. Among the list of pulsars matching with an eROSITA source, 52 have a counterpart in the 4FGL-DR4 catalog (see Sect. 3.3). For 15 investigated sources, we detected an X-ray counterpart for which we did not find a previous report in the literature. These sources are listed in Table 3 along with their basic pulsar parameters as adopted from the ATNF pulsar database (Manchester et al. 2005). Nine of these 15 pulsars have a γ-ray counterpart according to the 4FGL-DR4 catalog.

Of all 15 pulsars in Table 3, the X-ray counterpart of PSR J0837–2454 has the highest detection likelihood. This pulsar was discovered recently by Pol et al. (2021). Its age and distance estimates are 28.6 kyr and 0.2–0.9 kpc (Pol et al. 2021), respectively. Its X-ray counterpart is detected only in the 0.2–0.7 keV band, implying that the pulsar has a very soft X-ray spectrum. The 30 source counts detected in eRASS:4 do not support detailed spectral modeling due to the sparse photon statistics, but approximately constrain the interstellar column absorption toward PSR J0837–2454 to N_H ~ 10²⁰ cm⁻². This is in agreement with the column density through the Galaxy in the direction of PSR J0837–2454 which is N_H ≈ 7.5 × 10²⁰ cm⁻² (HI4PI Collaboration 2016). The low galactic column density further implies that the distance to PSR J0837–2454 is likely rather small. In any case, it is in stark contrast with the large dispersion measure toward the pulsar, which would imply a column density around N_H ≈ 4 × 10²¹ cm⁻² (Pol et al. 2021; He et al. 2013). This may imply an unusual ionization state of the interstellar medium along this line of sight, with an increased abundance of free electrons, meaning an enhanced ionization fraction.

Fig. 7 Example diagnostic plot for visual inspection. The top panels display eRASS:4 images of the position of the source 4FGL J1544.2−2554 in the energy bands 0.2–0.6, 0.6–2.3, and 2.3–5.0 keV, overlaid with the 1σ (blue dashed line) and 95% error ellipse of the 4FGL source (blue solid line), and the position of the detected X-ray source (red). The pixel size is 12″. The bottom panels show the local distribution of Gaia DR3 (left), CatWISE2020 (center), as well as SIMBAD entries and radio sources (right) around the X-ray source position (error bars). In the left and center panels, the size of the markers is proportional to the source brightness, and their color reflects the photometric colors (with bluer markers indicating bluer colors). In the left panel, pentagonal (triangular) markers correspond to sources with statistically significant (insignificant) proper motion or parallax, circular markers indicate the lack of a complete astrometric solution for a source.

5 Results and discussion

5.1 Catalog of candidate pulsar counterparts to unassociated 4FGL sources

Table 4 shows our 50 most likely pulsar candidates after visually inspecting all potential matches between eRASS:4 and 4FGL-DR4 catalogs. The full catalog of candidate matches, ranked by descending probability of pulsar-type matches $P_i^{\text{PSR}}$ is available in electronic form at CDS. This probability combines the probabilities of having a young and a recycled pulsar, meaning $P_i^{\text{PSR}}=P_i^{\text{YNG}}+P_i^{\text{MSP}}$. In addition, we also provide the visual inspection plots for each 4FGL source and a catalog of those sources that were removed from our final sample during visual inspections. Table B.1 gives an overview of the released catalog’s columns. Apart from the posterior match probabilities for different source classes, this catalog also includes X-ray properties such as source positions, errors, and fluxes (derived from count rates), as well as the individual quantities constructed to compute the match probabilities, such as flux ratios or γ-ray-based source classifications. Figure 9 summarizes the pulsar-match probabilities contained in our catalog and predicts how many new pulsars one may ideally expect to detect when observing a certain number of candidate counterparts, assuming our probabilities are quantitatively reliable (see Sect. 3.6). For instance, 6–11 pulsar detections may be expected among the top 20 of our match candidates or around 30–40 detections among the top 200. Among our candidates, an approximately equal number of MSPs and younger pulsars is predicted to be detected. While MSPs as high-latitude sources are usually in less crowded regions and therefore allow for more secure associations with their X-ray counterpart, unidentified sources at low latitudes, which are likelier to be young pulsars, are more numerous. This is substantiated by Fig. 10, which depicts the spatial distribution and match probabilities of potential pulsar-type counterparts to 4FGL sources over the western Galactic hemisphere.

To our knowledge, this work constitutes the first probabilistic cross-match between a Fermi-LAT source catalog and an X-ray catalog, and while only a few sources have pulsar-type match probabilities >50%, the expected yield of possible follow-up campaigns can be accurately quantified. Nonetheless, even though around 92% of all associated sources in 4FGL-DR4 are classified pulsars or blazars, the assumption that all unidentified sources belong to either of those source classes is most likely an oversimplification. This assumption neglects the existence of further γ-ray bright source types, such as supernova remnants, pulsar wind nebulae, or star-forming regions, which are particularly relevant for Galactic unassociated sources. Despite this limitation of our analysis, as Sect. 3.6 shows, our algorithm can efficiently identify pulsar-type counterparts to 4FGL sources, as verified on the set of associated sources. Hence, we believe our final candidate list offers a reliable starting point for any follow-up campaign.

Fig. 8 Verification of match probabilities. The top panels indicate the distribution of match probabilities $P_i^{\text{PSR}}/P_i^{\text{AGN}}$ (left/right) for associated 4FGL sources. The blue (orange) histograms indicate the match probabilities assigned to sources known to be the correct (incorrect) counterparts of the given 4FGL source, with the black dashed line indicating the sum of the two distributions. The lower panels show calibration curves for the two classes, plotting the observed fraction of correctly detected counterparts ƒdet, for sources binned by predicted match probability $P_i^{\text{PSR}}/P_i^{\text{AGN}}$. The red line indicates the 1 : 1 relation expected for an ideal, unbiased classifier.

Fig. 9 Overview of match probabilities. The left panel shows a histogram of pulsar−match probabilities $P_i^{\text{PSR}}$ (black, solid line) and the expected number of pulsar detections Nexp above a certain threshold of $P_i^{\text{PSR}}$ (black, dashed), after visual inspection. The green and blue dashed lines indicate the expected numbers of young pulsars and MSPs above the respective threshold. The right panel shows the dependence of the expected number of detections with the number of investigated X-ray sources Ninv for all pulsars and the two subsets in green and blue. The shaded regions indicate the expected statistical error on the predicted number.

5.2 Remarks on similar multiwavelength studies

Systematic counterpart searches for unidentified Fermi-LAT sources have been carried out with multiple instruments and at multiple wavelengths (see Smith et al. 2023). For instance, Bruzewski et al. (2023) identified a sample of high-likelihood pulsar candidates by searching for steep-spectrum radio sources coincident with 4FGL sources. Unfortunately, through a brief positional cross-match, we found that none of their sources were likely to be a counterpart to any of our X-ray sources.

Furthermore, a brief literature search yielded that, while a few of our X-ray sources have been known to be promising pulsar candidates for several years (Saz Parkinson et al. 2016; Dai et al. 2017), to our knowledge, none of our candidates have been independently confirmed to be pulsars yet. However, several likely candidates from a previous iteration of this project applied to the 4FGL-DR3 catalog (Abdollahi et al. 2022) have recently been confirmed to be very likely pulsars. This includes the redback system PSR J1803−6707 (Clark et al. 2023) or the MSP PSR J1259−8148 (Smith et al. 2023) discovered by the TRAPUM collaboration, but also several sources now classified as “binaries”, namely 4FGL J1702.7−5655, 4FGL J1910.7−5320, and 4FGL J1120.0–2204 (Corbet et al. 2022; Swihart et al. 2022a; Au et al. 2023). No pulsations have been detected yet for the latter sources, but periodic optical variability on hour timescales is observed, as is typical for MSP binaries.

Similar systematic studies to this work have been carried out using data from the Swift-XRT follow-up program of unasso-ciated Fermi sources (Stroh & Falcone 2013), for instance, the recent machine learning classification of X-ray counterparts to 4FGL sources by Kerby et al. (2021a). To compare the results of their method to ours, we positionally cross-matched their output catalog to our set of X-ray sources (before applying the pulsar-specific cuts in Sect. 3.4) with a radius of 15″. In Fig. 11, the output probabilities for matched sources are compared, considering the overall probability of a given match. Generally, secure matches of X-rays to γ-rays tend to agree between the two works, with all of these being classified as likely blazars. Similarly, all sources classified as more likely to be pulsars than blazars $\left(P_i^{\text{PSR}}>P_i^{\text{AGN}}\right)$ by Kerby et al. (2021a) are recovered quite well here. In two cases, however, sources deemed likely blazars in their work are classified as much more likely to be pulsars by us. In both cases (4FGL J1357.3−6123 and 4FGL J1639.8−4642c), the Fermi source lies in the Galactic plane, which in principle makes a pulsar scenario preferable. However, on the other hand, the X-ray fluxes of the sources observed by Swift are several orders of magnitude larger than those observed by eROSITA, and the source positions are 6″ and 13″ apart, respectively. Hence, it seems questionable whether truly identical sources were truly detected by Swift and eROSITA. If the vastly differing fluxes indeed stemmed from a variable single astrophysical source, this source would be unlikely to be a pulsar, as stable emission is expected on long timescales.

Importantly, the work of Kerby et al. (2021a) assumes that if an X-ray source is detected within an observed γ-ray error ellipse, it is in fact its likely counterpart. While this approach appears suitable for X-ray-bright sources such as blazars, in the case of γ-ray pulsars, it is important to consider the possibility of not having detected the correct counterpart at all (see Fig. 11), as given the expected X-ray fluxes, in most cases an X-ray non-detection is a priori more likely than a detection (see Sect. 3.3). Quantitatively estimating the match probability requires knowledge of the density and properties of the background X-ray source population, which underlines why an all-sky survey is an ideal tool to carry out such a task.

Fig. 10 Distribution of the best pulsar-type match candidates to Fermi sources in the western Galactic hemisphere. The color bar indicates the maximum probability of having an X-ray-detected pulsar counterpart $P_i^{\text{PSR}}$ for each source. Sources deemed likelier to be young pulsars are marked with squares and likelier MSPs with triangles. Sources that did not pass visual inspection are marked with circles.

Fig. 11 Comparison between our results and those of Kerby et al. (2021a) for the sample of X-ray sources common to both works. We plot the relative probability for a source to be a pulsar rather than a blazar, as determined here against the results from Kerby et al. (2021a). The marker colors indicate the total probability of an X-ray source being the correct counterpart to the matched Fermi source $P_i=P_i^{\text{PSR}}+P_i^{\text{AGN}}$.

5.3 Prospects for multiwavelength follow-up

The results of our cross-matching effort do not constitute immediate identifications of pulsars, as the positional uncertainties of the 4FGL γ-ray sources are much too large to perform secure associations. However, we believe that our results are ideally suited as input to multiwavelength follow-up campaigns with the target of detecting new pulsars (see Smith et al. 2023), especially since a quantification of the expected success rate is possible with our output probabilities.

The most natural choice for such a follow-up campaign would be a radio search for pulsed emission, similar to Clark et al. (2023). In particular, our X-ray positions, accurate to around 5″, would vastly reduce the area, or the number of barycentric corrections, to be searched compared to the arcminute-sized error ellipses of 4FGL sources. Furthermore, a search for pulsar candidates in the imaging domain via deep pointed observations similar to Bruzewski et al. (2023) could be a viable alternative. While this would not allow for directly detecting pulsations, the combination of X-ray emission and a steep-spectrum radio source coincident with a γ-ray source would make a rotation-powered pulsar appear very likely.

For those X-ray sources that we classify as high-likelihood MSP candidates, optical observations with time-resolved photometry would likely constitute a similarly effective path, which may be somewhat less expensive than radio follow-up. Since, due to their evolutionary history, MSPs are expected to primarily exist in binary systems, optical emission is expected from their companions. Due to their short orbital periods, in particular for the spider subclass (Roberts 2013), modulations of the optical emission are expected on time scales of hours. While optical data do not allow for the detection of pulsations from the pulsar magnetosphere itself, the detection of such rapid periodic modulation from an X-ray bright source constitutes a “smoking gun” for a binary pulsar and is an ideal tool for the identification of (candidate) MSP binaries (e.g., Clark et al. 2023; Corbet et al. 2022; Swihart et al. 2022b,a; Romani et al. 2012; Romani & Shaw 2011).

Finally, while the majority of MSPs is radio-loud, around half of the young γ-bright pulsars are radio-quiet (Smith et al. 2023). For this class of objects, neither optical nor radio searches would appear particularly promising. Hence, one possible, yet observationally expensive, avenue would be to carry out X-ray follow-up campaigns with the target of directly detecting pulsations. Given the X-ray faint nature of our typical pulsar candidates in this work, XMM-Newton (Strüder et al. 2001; Turner et al. 2001) appears likely to be the only current X-ray mission feasible for this purpose, since a high number of source counts is as crucial for detecting pulsations as the ability to separate source from background emission, which makes the Chandra (Weisskopf et al. 2002) or NICER (Gendreau et al. 2016) telescopes appear less attractive. However, even with XMM-Newton, this effort would require quite long exposures. For instance, a typical eRASS:4 source flux of F_X = 2 ×10⁻¹⁴ erg s⁻¹ cm⁻², corresponding to an XMM 0.2–12 keV count rate of 1.8 × 10⁻²ct s⁻¹ in a typical aperture, together with a pulsed fraction of 20%, would require an exposure of around 35 ks for a 5σ detection of a pulsed signal. The only viable alternative to X-ray follow-up would be to directly search the available γ-ray data from Fermi-LAT for pulsations at the candidate positions. Blind searches for γ-ray pulsations from unidentified Fermi-LAT sources were carried out successfully in the past via the distributed computing system Einstein@Home (Clark et al. 2017). These searches tend to suffer from having to search multiple possible positions due to the localization inaccuracy of Fermi-LAT, a problem which would be alleviated by a targeted search of candidate positions with arcsecond-level accuracy (such as Smith et al. 2019).

6 Summary and conclusion

In this work, we have conducted a probabilistic cross-matching analysis between the 14-yr catalog of Fermi-LAT γ-ray sources (4FGL-DR4), and the SRG/eROSITA X-ray catalog based on four all-sky surveys. The primary objective of this effort was the identification of new high-energy rotation-powered pulsars by searching for likely X-ray counterparts to unassociated γ-ray sources. This was motivated by the fact that pulsars are the second most common class of identified 4FGL sources behind blazars. A dedicated multiwavelength study of the eROSITA-detected population of (candidate) blazars will be published elsewhere (S. Hämmerich et al., in prep.). Given the large positional uncertainty inherent to Fermi-LAT sources, the identification of an X-ray counterpart facilitates their identification, exploiting the characteristic thermal surface or magnetospheric X-ray emission from pulsars, even if little more than arcsecond-level positional information is provided.

For the identification of new pulsar candidates, we employed a Bayesian cross-matching scheme designed for astronomical catalogs, which combines positional overlap with prior information on the likelihood of a given match to produce a probabilistic statement on the quality of a given candidate match of pulsar- or blazar-type. More specifically, we constructed a random forest classifier, which, based on the γ-ray properties of the investigated Fermi-LAT sources, predicts the probability for it to be of a blazar or pulsar (young or recycled) nature. While this classifier achieves similar accuracies as in previously published works, it predicts a much larger fraction of pulsars among unassociated sources. We attribute this to our emphasis on a balanced treatment of the source classes and the exclusion of explicit source fluxes from the approach.

These ingredients were combined with a data-driven constraint on the expected X-ray flux for a match of either type based on the observed γ- to X-ray flux ratios for sources with known multiwavelength counterparts. This effort clearly confirmed previously established results, namely the higher relative X-ray brightness of blazars, and the large spread in X-ray fluxes of young pulsars. Finally, the X-ray images of our resulting candidates were visually screened in order to identify potential spurious detections. In addition, we investigated optical and infrared catalogs with the target of excluding sources classifiable as non-pulsar based on their multiwavelength counterparts.

We verified our method of probabilistic cross-matching by applying it “blindly” to the set of previously associated sources in 4FGL-DR4, finding decent agreement between the predicted match probabilities and the fraction of counterparts that were in fact correctly identified. Our final list of candidate matches contains around 900 X-ray sources with probabilities of being the X-ray counterpart to a γ-ray pulsar hidden in an unassociated 4FGL source $P_i^{\text{PSR}}\ge 0.02$. Around 50 sources exhibit probabilities $P_i^{\text{PSR}}\ge 0.02$, most of which lie close to the Galactic plane. Based on our posterior match probabilities, we predict that a “perfect pulsar detector” should find 30–40 new pulsars among our 200 best candidates, with about equal shares of young energetic and old millisecond pulsars. Support for the reliability of our method comes from several sources now known to be pulsars, which were likely candidates in previous iterations of this project, and from its qualitative agreement with X-ray follow-up campaigns by Swiƒt-XRT.

In addition to our counterpart search for Fermi-LAT sources, using our homogeneous survey coverage of the X-ray sky, we have performed a systematic search for X-ray counterparts of known rotation-powered pulsars, including all known radio and γ-ray pulsars. Thereby, we have identified 15 previously unknown X-ray counterparts of pulsars, including the soft emission from the recently discovered young and nearby PSR J0837–3754.

We hope that our study will motivate new, or support existing, multiwavelength identification programs of pulsars in unassociated Fermi-LAT sources. Targeted searches for pulsations in the radio or γ-ray domain could profit from our candidate positions with an arcsecond-level accuracy, as it would imply a large reduction in the area to search, relative to the γ-ray error ellipses. Similarly, optical searches for periodically variable counterparts, induced by the orbital modulation of the brightness of a potential binary companion to the pulsar, would be facilitated by only having to search the X-ray uncertainty region rather than the full γ-ray error ellipse, which may contain thousands of sources.

Acknowledgements

We would like to thank Dusán Tubín-Arenas for providing us with X-ray upper limits for all associated 4FGL sources. Further, we would like to thank Julien Wolf and Mara Salvato for useful discussions on multiwavelength counterparts in early phases of this project. W.B. acknowledges fruitful discussions with Luciano Nicastro. M.G.F.M. acknowledges support by the International Max-Planck Research School on Astrophysics at the Ludwig-Maximilians University (IMPRS). This work is based on data from eROSITA, the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tubingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA. The eROSITA data shown here were processed using the eSASS software system developed by the German eROSITA consortium. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular, the institutions participating in the Gaia Multilateral Agreement. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018). Further, we acknowledge the use of the Python packages Matplotlib (Hunter 2007), SciPy (Virtanen et al. 2020), and NumPy (Harris et al. 2020).

Appendix A Tuning of random forest classifiers

The goal of classifying the 4FGL-DR4 catalog into likely blazars and pulsars based on machine learning suffers from a strong inherent imbalance: since about ten times more identified blazars can be used for training the classifier, a standard random forest would tend to learn a bias against the minority class, which happens to be the class of interest, here. Hence, in order to assign equal “prior” probabilities to the AGN and pulsar scenarios, samples of similar size need to be used for training. To achieve this, one could, on one hand, try to oversample the minority class, either by using sources multiple times or by producing artificial samples with the synthetic minority oversampling technique (SMOTE; Chawla et al. 2002). A somewhat more natural choice is to undersample the majority class, meaning to draw only a random subsample for each individual decision tree, so that equal sample sizes are effectively given (Chen et al. 2004). As we found that, even for low recall, the SMOTE method struggled to achieve a pure pulsar sample in our tests, we decided to employ the latter method, which is implemented in the BalancedRandomForest classifier in the imbalanced–learn Python package (Lemaître et al. 2017).

Generally, we found that the hyperparameters of our classifier had a smaller impact on its performance than the choice and definition of the input features. The classifiers discussed in the following employed a maximum tree depth of ten and a maximum number of features per node of two, which was one of the setups found to yield optimal performance. In order to allow for reliable probabilistic statements in our output, we trained a large number of trees, that is, 10 000.

Fig. A.1 Distribution of 4FGL-DR4 sources in feature space. In each panel, we show the distribution of pulsars (blue), blazars (red), and unassociated sources (black) in different parameters. The left panel compares the observed energy flux Fγ in the 100 MeV – 100 GeV band to its estimated error σγ. The central panel compares the significance of detection with that of spectral curvature for a log-parabola model. In the right panel, this is compared with the measure of variability on monthly to yearly scales. In each panel, the dashed gray line indicates the direction perpendicular to which we defined our composite parameters to be able to better separate the source types.

Fig. A.2 Illustration of recursive feature elimination. The top panel illustrates the evolution of the average precision of our classifier with the number of features. The bottom panel illustrates the evolution of the importance of individual features, as the least important feature is recursively removed. The shaded bar indicates the set of features we found to be optimal. The curves indicate the parameters numbered by their reverse order of elimination: (1) SigCombined, (2) ModSigCurv, (3) LP_beta, (4) LP_SigCurv, (5) K24, (6) SymLat, (7) Frac_Variability, (8) EFluxErr, (9) HR34, (10) Variability_Index, (11) HR45, (12) PL_Index, (13) HR 12, (14) DeltaInd, (15) SymLon, (16) K13, (17) Pivot_Energy, (18) LP_Index, (19) K35, (20) HR23.

Several possible biases may impact the performance of our binary classifier, for instance, a possible covariate shift, where the distribution of input features is different between training and target samples (Malyshev 2023). A natural example of this is the source flux, since the sample of unassociated sources is systematically fainter in γ-rays than the set of associated sources, in particular pulsars (see Fig. A.1). This is because the direct detection of γ-ray pulsations requires a comparatively large number of photons, hence a sufficiently bright γ-ray source (Saz Parkinson et al. 2016). In order to counteract biases in our classification against the pulsar class introduced by this effect, we did not include any absolute fluxes in our set of candidate features. Further, we constructed composite features, in which we attempted to remove the dependence on source flux or significance.

The three panels in Fig. A.1 illustrate this interdependence between some features which have been shown to be important for classification (e.g., Salvetti 2014; Saz Parkinson et al. 2016). Clearly, the flux error strongly scales with the measured energy flux, even though the lower envelope exhibits a slightly shallower scaling than expected in a purely Poissonian regime with around ${\sigma }_{\gamma }\propto F_{\gamma }^{0.4}$. The scatter relative to this relation is probably related to varying background levels in the vicinity of the sources. Similarly, the curvature significance for the log-parabola model LP_SigCurv scales with the source detection significance, but pulsars clearly tend to show more significant curvature in a relative sense (Abdollahi et al. 2022). Similarly, due to their dependence on source flux, LP_SigCurv exhibits a weak scaling with the measure of variability, most clearly for bright blazars. We constructed several composite features based on the observed scalings. Outside the modified flux error measure (see left panel in Fig. A.1), the primary goal of these was to optimally separate the two source classes in the displayed parameter spaces, rather than being strictly flux-independent. Furthermore, similarly to previous works (Bhat & Malyshev 2022; Germani et al. 2021; Saz Parkinson et al. 2016, e.g.,), we constructed several empirical hardness ratios as well curvature measures based from the band fluxes given in the catalog. Our full list of candidate features is given in Table A.1.

We followed the established approach of recursive feature elimination (Luo et al. 2020) for identifying the optimal set of features to use for our classification task: Starting from a classifier with all candidate features, we repeatedly determined its performance, characterized by its average precision measured with repeated five-fold cross-validation, and successively removed the feature with the lowest importance, as illustrated in Fig. A.2. The average precision measures the ability of the classifier to produce a pure and complete sample of pulsars in the presence of a class imbalance, hence we chose the version of our classifier with 12 features (see Table 1), which maximizes this quantity. Many of the features with the highest importance for the classifier are directly related to the spectral curvature, and in a few cases also to the relative variability of the sources. This is an unsurprising behavior, as it has been known for a while that pulsars show less variability and have spectra which are more curved than blazars (Salvetti 2014; Saz Parkinson et al. 2016).

Fig. A.3 Receiver operating characteristic of the binary classifier separating young pulsars and MSPs. This graph compares the true positive rate, meaning the fraction of young pulsars that are correctly identified, to the false positive rate, meaning the fraction of MSPs wrongly identified as young pulsars, dependent on the decision threshold.

As a final step, we converted the prediction of our random forest into a quantitative probability $P_{\gamma }^{\text{PSR}}$, specifying the chance for a given γ-ray source to be a pulsar, assuming both scenarios are a priori equally likely. To do this, we used the output probabilities of our cross-validation and created calibration curves. This means we computed the fraction of pulsars f_PSR in the output sample dependent on the probability p assigned by the classifier. Since the number of blazars N_AGN is much larger than the number of pulsars N_PSR in the test sample, corrected this value to obtain the “balanced” fractions which would be observed for equally abundant classes: $$\frac{f_{\text{PSR}}^{\text{bal}}}{1-f_{\text{PSR}}^{\text{bal}}}=\frac{N_{\text{AGN}}}{N_{\text{PSR}}}\frac{f_{\text{PSR}}}{1-f_{\text{PSR}}}.$$(A.1)

Given the binary nature of our classifier, a logistic curve of the form $$f(x\mid a,b)=\frac{1}{1+\exp (-ax+b)},$$(A.2)

where x is related to the classifier output p as x = ln p/ (1 − p), is a natural choice for its calibration. By fitting this function to the observed dependence of pulsar fraction $f_{\text{PSR}}^{\text{bal}}$ on classifier prediction p, we obtained a quantitatively calibrated estimator $P_{\gamma }^{\text{PSR}}=f(p)$ for the pulsar probability of any given γ-ray source.

Analogously to our classification of pulsars versus blazars, we created a binary classifier for estimating the probability for a given source to be a young pulsar, rather than a millisecond pulsar, assuming its pulsar nature. We followed the same approach as for our primary classifier, with the only exception that we used the receiver operating characteristic (ROC), to measure the performance of our classifier, as neither class is of larger interest for our study and source fractions are about balanced, here. The optimal ROC curve, with an area under the curve of 0.964, is displayed in Fig. A.3. The optimal set of features for this classifier are shown in Table A.2, and are fundamentally different from the pulsar-blazar case. The highest importance was ascribed to the Galactic latitude of the source, as young pulsars tend to be much closer to their birth site, meaning in the Galactic plane. Further important features include the relative flux errors, which is probably related to the crowdedness of the field, and the overall slope of the γ-ray spectrum, as MSPs typically have somewhat harder spectra (Abdo et al. 2013).

Finally, we note that we also tested whether the two classifiers used in our work could be combined into a single three-way classifier separating blazar-like sources, young pulsars, and MSPs. However, we found that the performance of this single classifier fell slightly short of that combining the two individual ones. Quantitatively, the three-way classifier achieved an average precision of 0.916 in the separation of pulsars from AGN, and a ROC AUC of 0.959 in the subclassification of pulsars into young and recycled pulsars, compared to scores of 0.920 and 0.964 by the two separate classifiers, respectively.

Appendix B Supplementary tables

Here, we provide additional information for readers who wish to reproduce the analysis performed in this work, or employ our results for their own research. Table B.1 describes the columns of our full catalog of candidate pulsar counterparts to 4FGL sources, which is provided in electronic form at CDS. Table B.2 lists all X-ray detections of pulsars in the 4FGL catalog, along with the X-ray fluxes derived from their count rates, as displayed in Fig. 4.

References

All Tables

All Figures

	Fig. 1 Precision-recall curve for the binary classifier separating pulsar-and blazar-candidate γ-ray sources. This graph compares the recall of the classifier, i.e., the fraction of pulsars in the test sample that are correctly identified, to its precision, i.e., the fraction of true pulsars among all objects selected as pulsars, dependent on the decision threshold.
In the text

Fig. 2 Distribution of binary classification probabilities $P_{\gamma }^{\text{PSR}}$ for the targeted unassociated 4FGL sources. This is compared to the probabilities evaluated on PSR- and AGN-type “out-of-box” test samples during cross-validation, meaning targets the classifier did not see during training. The fractions of misclassified AGN and pulsars in the test sample are 4.2% and 4.9%, respectively.

In the text

	Fig. 3 Distribution of binary classification probabilities $P_{\gamma }^{\text{PSR}}$, indicated by the color scale of the markers, of unassociated 4FGL sources across the Galaxy. The gray shading marks the eastern Galactic hemisphere, for which no eROSITA X-ray data are available.
In the text

Fig. 4 Distribution of X-ray (0.2–2.3 keV) and γ-ray (100 MeV – 100 GeV) fluxes of young pulsars (blue circles) and MSPs (green). One-sigma X-ray upper limits are indicated with triangles in cyan (young pulsars), lime green (MSPs), and magenta (AGN), respectively. In order to provide a clean plot, the distribution of detected AGN is illustrated by the red filled contours, marking its one, two, and three-sigma outlines.

In the text

	Fig. 5 Distribution of the γ-to-X-ray flux ratio of 4FGL young pulsars (top), MSPs (middle), and AGN (bottom). In both panels, the distribution of the flux ratio of detected sources (dotted histogram) and corresponding one-sigma lower limits (dashed histogram) are indicated, together with a best-fitting log-normal distribution (solid line).
In the text

Fig. 6 Example for the computation of the expected X-ray flux for a specific γ-ray source, 4FGL J0159.8–2234. The top panel indicates the normalized X-ray flux distributions of field sources (black), and for the assumption of a young pulsar, MSP, or AGN nature. The bottom panel depicts the corresponding ratio of probabilities. The markers indicate the fluxes of the detected X-ray sources around the Fermi source. For reference, the γ-ray flux of the source is log Fγ = −12.00, and the limiting X-ray flux is λ = −13.95.

In the text

Fig. 7 Example diagnostic plot for visual inspection. The top panels display eRASS:4 images of the position of the source 4FGL J1544.2−2554 in the energy bands 0.2–0.6, 0.6–2.3, and 2.3–5.0 keV, overlaid with the 1σ (blue dashed line) and 95% error ellipse of the 4FGL source (blue solid line), and the position of the detected X-ray source (red). The pixel size is 12″. The bottom panels show the local distribution of Gaia DR3 (left), CatWISE2020 (center), as well as SIMBAD entries and radio sources (right) around the X-ray source position (error bars). In the left and center panels, the size of the markers is proportional to the source brightness, and their color reflects the photometric colors (with bluer markers indicating bluer colors). In the left panel, pentagonal (triangular) markers correspond to sources with statistically significant (insignificant) proper motion or parallax, circular markers indicate the lack of a complete astrometric solution for a source.

In the text

Fig. 8 Verification of match probabilities. The top panels indicate the distribution of match probabilities $P_i^{\text{PSR}}/P_i^{\text{AGN}}$ (left/right) for associated 4FGL sources. The blue (orange) histograms indicate the match probabilities assigned to sources known to be the correct (incorrect) counterparts of the given 4FGL source, with the black dashed line indicating the sum of the two distributions. The lower panels show calibration curves for the two classes, plotting the observed fraction of correctly detected counterparts ƒdet, for sources binned by predicted match probability $P_i^{\text{PSR}}/P_i^{\text{AGN}}$. The red line indicates the 1 : 1 relation expected for an ideal, unbiased classifier.

In the text

Fig. 9 Overview of match probabilities. The left panel shows a histogram of pulsar−match probabilities $P_i^{\text{PSR}}$ (black, solid line) and the expected number of pulsar detections Nexp above a certain threshold of $P_i^{\text{PSR}}$ (black, dashed), after visual inspection. The green and blue dashed lines indicate the expected numbers of young pulsars and MSPs above the respective threshold. The right panel shows the dependence of the expected number of detections with the number of investigated X-ray sources Ninv for all pulsars and the two subsets in green and blue. The shaded regions indicate the expected statistical error on the predicted number.

In the text

Fig. 10 Distribution of the best pulsar-type match candidates to Fermi sources in the western Galactic hemisphere. The color bar indicates the maximum probability of having an X-ray-detected pulsar counterpart $P_i^{\text{PSR}}$ for each source. Sources deemed likelier to be young pulsars are marked with squares and likelier MSPs with triangles. Sources that did not pass visual inspection are marked with circles.

In the text

Fig. 11 Comparison between our results and those of Kerby et al. (2021a) for the sample of X-ray sources common to both works. We plot the relative probability for a source to be a pulsar rather than a blazar, as determined here against the results from Kerby et al. (2021a). The marker colors indicate the total probability of an X-ray source being the correct counterpart to the matched Fermi source $P_i=P_i^{\text{PSR}}+P_i^{\text{AGN}}$.

In the text

Fig. A.1 Distribution of 4FGL-DR4 sources in feature space. In each panel, we show the distribution of pulsars (blue), blazars (red), and unassociated sources (black) in different parameters. The left panel compares the observed energy flux Fγ in the 100 MeV – 100 GeV band to its estimated error σγ. The central panel compares the significance of detection with that of spectral curvature for a log-parabola model. In the right panel, this is compared with the measure of variability on monthly to yearly scales. In each panel, the dashed gray line indicates the direction perpendicular to which we defined our composite parameters to be able to better separate the source types.

In the text

Fig. A.2 Illustration of recursive feature elimination. The top panel illustrates the evolution of the average precision of our classifier with the number of features. The bottom panel illustrates the evolution of the importance of individual features, as the least important feature is recursively removed. The shaded bar indicates the set of features we found to be optimal. The curves indicate the parameters numbered by their reverse order of elimination: (1) SigCombined, (2) ModSigCurv, (3) LP_beta, (4) LP_SigCurv, (5) K24, (6) SymLat, (7) Frac_Variability, (8) EFluxErr, (9) HR34, (10) Variability_Index, (11) HR45, (12) PL_Index, (13) HR 12, (14) DeltaInd, (15) SymLon, (16) K13, (17) Pivot_Energy, (18) LP_Index, (19) K35, (20) HR23.

In the text

	Fig. A.3 Receiver operating characteristic of the binary classifier separating young pulsars and MSPs. This graph compares the true positive rate, meaning the fraction of young pulsars that are correctly identified, to the false positive rate, meaning the fraction of MSPs wrongly identified as young pulsars, dependent on the decision threshold.
In the text