Searching for X-ray counterparts of unassociated Fermi -LAT sources and rotation-powered pulsars with SRG /eROSITA

,


Introduction
The Fermi Large Area Telescope (LAT) is by far the most sensitive γ-ray telescope in the GeV 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 years of data and supersedes three previous iterations of the fourth Fermi-LAT source catalog (Abdollahi et al. 2020;Ballet et al. 2020;Abdollahi et al. 2022).It contains 7195 γ-ray emitters from 50 MeV to 1 TeV, distributed over the whole sky.The two most prominent contributing source classes are blazartype 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, does not yet have a convincing association at lower energies, providing considerable discovery potential in the ⋆ e-mail: mgf.mayer@fau.deFermi-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 mismodelled 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 is found in binary systems with an optical companion, due to their formation requir-ing the 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 X / Ė ∼ 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 GeV γ-ray efficiencies are significantly larger, at L X / Ė ∼ 10 −2 -1 (Abdo et al. 2013).For Fermi-detected pulsars in particular, the GeV γ-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(Marelli et al. , 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 rotationpowered pulsars are expected to show detectable X-ray emission.The positions of X-ray sources are usually accurate on an arcsecond-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) 1 .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 observatory (SRG; 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 years 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.3keV 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 data set (named eRASS:4) reaching a typical point source sensitivity of around 1 − 2 × 10 −14 erg s −1 cm −2 in the main band.This paper employs the eRASS:4 data of the western Galactic hemisphere, i.e. in the range 180 • < l < 360 • .This data set 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 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 Xray 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 followup prospects.Finally, we summarize our work in Sect.6.

Input Data from Fermi-LAT and eROSITA
The target of our analysis is the set of all unidentified sources in the Fermi-LAT 14-year source catalog (4FGL-DR4, in the following 4FGL), which contains 7195 GeV γ-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. 2 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.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.

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 −2 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 multiwavelength counterparts (Sect.3.5).Finally, we tested our method using the known associations of 4FGL sources (Sect.3.6).

Classification of γ-ray sources
It has been shown numerous times (Salvetti 2014;Saz Parkinson et al. 2016;Kaur et al. 2019;Kerby et al. 2021a,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-learningbased classification of unassociated 4FGL-DR4 sources, with the main target of estimating a prior probability P 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 imbalanced nature of the training set (with a factor of ten more blazars than pulsars), we decided to use a BalancedRandomForest classifier (Chen & Breiman 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 AUC3 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 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 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).Fig. 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 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 5fold 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 PSR γ and P YNG γ , as defined in this section.
Article number, page 4 of 21 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 light red (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.

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 where ν 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 σ 1 and σ 2 is given by (see eqs. 16 & 28 in Budavári & Szalay 2008): (2) Hence, following from Bayes' theorem, the posterior probability of a given match is proportional to (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 locationdependent one ν = ν(x), which we estimated locally within a two-degree radius of each Fermi source.This yields the follow-ing expression for the astrometric component of our posterior: Here, we have decomposed the separation vector ψ into its x and y components.The positional covariance matrix components for 4FGL sources (ρ 1 , σ 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 σ 2 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).

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(Marelli et al. , 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_14 ), and γ-ray flux in the 100 MeV − 100 GeV band.In case of no matching source, we used the internal eROSITA upper limit server (see Tubín-Arenas et al. 2023) to determine the counts and background level within a 30 ′′ aperture5 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 the appendix 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 ∼ N(µ, σ).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 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. 2023) 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 fitted a simple empirical model for the background flux distribution ϕ B of the type: 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 classdependent 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 : 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 −12 erg s −1 cm −2 ).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 Xray 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.

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: where the index i runs over all possibly associated X-ray sources, and t = {YNG, MSP, AGN} indicates the considered object types.Here, ρ t i can be understood as an unnormalized match probability, taking into account the γ-ray properties (via the source-type prior P 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 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 t any = i P t i , and of having any match of any type by P any = i t P t i .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.
After having computed the "pulsar-match" probability 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 PSR i ≥ 0.02.This left us with around 1100 candidate matches, 200 of which have a probability P PSR i ≥ 0.1.

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 multiwavelength 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 et al. 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 frac-tion 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 optically fainter 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 6 .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.

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 PSR γ , using our out-of-bag estimates created via 5-fold crossvalidation (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 PSR i or P AGN i , and evaluating the number of correctly and incorrectly identified counterparts in each bin, we 6 https://zenodo.org/doi/10.5281/zenodo.10634392are 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 PSR i /P AGN i , and the observed fraction of correctly identified counterparts f det , despite a larger scatter for the pulsar class due to the much smaller test sample.

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:1 showed that the astrometric accuracy varies between ∼ 4 − 15 arcsec 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 data set, we found that 97 Notes.The full list of sources is available online at CDS.The columns α X , δ X , σ X give the position and uncertainty of the matched X-ray source; L det and L ext specify its detection and extent likelihood (Brunner et al. 2022).P PSR γ gives the prior probability of the respective γ-ray source to correspond to a pulsar (Sect.3.1).Finally, P PSR i (P AGN i ) indicates the Bayesian probability for a particular X-ray source to be the correct match and for the source's nature to be a pulsar (blazar), respectively.P YNG 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 in-vestigated 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.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 towards PSR J0837−2454 to N H ∼ 10 20 cm −2 .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 20 cm −2 (HI4PI Collaboration et al. 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 strong contrast with the large dispersion measure toward the pulsar, which would imply a column density around N H ≈ 4 × 10 21 cm −2 (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.

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 PSR i is available in electronic form at CDS.This probability combines the probabilities of having a young and a recycled pulsar, meaning 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 γ-raybased source classifications.Fig. 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 followup 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 followup campaign.

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 Xray 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 unassociated 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 pulsarspecific 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 (P PSR i > P AGN i ) 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 nondetection 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.

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 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 exist in binary systems primarily, optical emission is expected from their companions.Due to their short orbital periods, particularly 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 are 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 −14 erg s −1 cm −2 , corresponding to an XMM 0.2−12 keV count rate of 1.8×10 −2 ct s −1 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).

Summary and conclusion
In this work, we have conducted a probabilistic cross-matching analysis between the 14-year catalog of Fermi-LAT γ-ray sources (4FGL-DR4), and the SRG/eROSITA X-ray catalog based on four all-sky surveys.The primary target 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 eROSITAdetected 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 blazar or (young or recycled) pulsar 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 nonpulsar 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 PSR i ≥ 0.02.Around 50 sources exhibit probabilities P PSR i ≥ 0.2, 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 Swift-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 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.

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 & Breiman 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 10 and a maximum number of features per node of 2, 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, i.e. 10 000.
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 σ γ ∝ F 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 logparabola 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).
As a final step, we converted the prediction of our random forest into a quantitative probability P 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: 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 bal PSR on classifier prediction p, we obtained a quantitatively calibrated estimator P 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).

Fig. 2 :
Fig. 2: Distribution of binary classification probabilities P 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.
See Appendix A for the reasoning behind the definition of the compound features.(a)Catalog column (seeAbdollahi et al. 2022)

Fig. 3 :
Fig.3: Distribution of binary classification probabilities P PSR γ , indicated by the color scale of the markers, of unassociated 4FGL sources across the Galaxy.The grey shading marks the eastern Galactic hemisphere, for which no eROSITA X-ray data are available.

Fig. 5 :
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 :
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.

Fig. 7 :
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.

Fig. 8 :
Fig. 8: Verification of match probabilities.The top panels indicate the distribution of match probabilities P PSR i /P AGN i (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 f det , for sources binned by predicted match probability P PSR i /P AGN i .The red line indicates the 1 : 1 relation expected for an ideal, unbiased classifier.
for the young pulsar and MSP subtypes.The entries have been sorted by descending P PSR i . (a) Identified as likely MSP in Dai et al. (2017); (b) Possible stellar counterpart; (c) Star-rich region; (d) Possible binary counterpart; (e) Classified as likely blazar by Kerby et al. (2021a); (f) Possible Be star counterpart; (g) Classified as likely pulsar by Kerby et al. (2021a) .

Fig. 9 :
Fig. 9: Overview of match probabilities.The left panel shows a histogram of pulsar-match probabilities P PSR i (black, solid line) and the expected number of pulsar detections N exp above a certain threshold of P PSR i (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 N inv for all pulsars and the two subsets in green and blue.The shaded regions indicate the expected statistical error on the predicted number.

Fig. 10 :
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 PSR i 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 :
Fig.11: Comparison between our results and those ofKerby 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 fromKerby 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 PSR i + P AGN i .

Fig
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.

Table 1 :
Optimal set of features and associated importance for the classifier separating pulsar-like and blazar-like γ-ray sources.

Table 2 :
Results of fitting flux ratios F γ /F X with a log-normal distribution N(µ t , σ t ) for different source classes t.Notes.Errors are given at a one-sigma level.

Table 3 :
Observational properties of the ATNF pulsars with newly discovered X-ray counterparts in eRASS:4 data.
(Brunner et al. 2022005)umns are as follows: PSR: Pulsar name, γ-ray counterpart: Name of the source in the 4FGL catalog.The pulsar parameters P (rotation period), Ṗ (period derivative), d (distance), (spin-down) age, and Ė (spin-down power) are taken from the ATNF pulsar database(Manchester et al. 2005).Rate is the count rate in the 0.2 − 2.3 keV band.L det is the detection likelihood for the 0.2 − 2.3 keV band(Brunner et al. 2022).The list is sorted in descending order of Ė/4πd 2 .

Table 4 :
List of top 50 candidate pulsar-type matches to 4FGL sources.

Table A .
2: Optimal set of features and their importance for the classifier separating young pulsars and MSPs among γ-ray sources.