Candidate isolated neutron stars in the 4XMM-DR10 catalogue of X-ray sources

ABSTRACT

1 INTRODUCTION

Neutron stars are the remnants of massive stars whose cores collapse during the supernova explosions. They are usually detected as radio (∼2500 known, Manchester et al. 2005) and/or gamma-ray (∼300, Abdollahi et al. 2020) pulsars, thanks to their beamed non-thermal emission. In addition, there is a group of isolated neutron stars (INSs) from which pulsed non-thermal emission is not detected.

This could be caused by an unfavourable orientation of the rotation and magnetic axes, or by the intrinsic lack, or faintness, of non-thermal magnetospheric emission. On the other hand, these objects can be discovered through their thermal X-ray emission, which arises from the cooling neutron star surface and, being emitted nearly isotropically, can be seen independently of the star orientation. This thermal emission is characterized by a soft X-ray spectrum that peaks between 0.2 and 2 keV (and references therein Potekhin et al. 2020) and gives only a small contribution in the optical band, leading to high values of the X-ray to optical flux ratio, F_X/F_O ≳ 10³. Pulsations with broad, nearly sinusoidal profiles, and usually small pulsed fractions (e.g. Tiengo & Mereghetti 2007), can be detected. These are likely caused by a non-uniform surface temperature and/or beaming effects due to the presence of a magnetized atmosphere.

The prototype of these thermally emitting INSs was discovered with the ROSAT satellite (Walter, Wolk & Neuhäuser 1996), and, subsequently, six other similar sources were found in the ROSAT All-Sky Survey (RASS, Voges et al. 1999). These sources, generally referred to as X-ray-dim isolated neutron stars (XDINSs), have distances in the range of 100–500 pc, X-ray spectra with blackbody temperatures of 45–110 eV, and spin periods of a few seconds (Haberl 2007; Turolla 2009).

Thanks to their high sensitivity, X-ray satellites, such as XMM–Newton and Chandra, should be able to detect dimmer (and thus farther) thermally emitting INSs. However, they provide a smaller sky coverage than all-sky surveys. Furthermore, absorption by the interstellar medium strongly suppresses soft X-rays, thus reducing the possibility of detecting farther and more absorbed objects. In fact, further searches using these satellites led only to a small increase in this sample. For example, Pires, Motch & Janot-Pacheco (2009b) analysed the 2XMMp catalogue of serendipitous sources discovered with XMM–Newton and found a few possible INS candidates, among which the most promising is 2XMM J104608.7−594306 (Pires et al. 2009a, 2015).

Based on a population synthesis model accounting for the distribution of neutron star birthplaces and of the local (<3 kpc) interstellar absorption, Posselt et al. (2008) suggested that the most promising regions to look for INSs are in the direction of rich OB associations. The eROSITA/SRG instrument (Predehl et al. 2021) has recently completed the first of the planned 4 yr of its all sky survey. This has provided the first coverage of the whole sky carried out in the soft X-rays after the RASS, although in a slightly different energy range: 0.1–10 keV with respect to 0.1–2.4 keV. Predictions for the expected number of INS in the complete eROSITAsurvey have been reported by Pires, Schwope & Motch (2017), who estimate about 90 discoveries, down to 0.2–2 keV fluxes of ∼10⁻¹⁴ erg s⁻¹ cm⁻² in the whole sky. However, their secure identification in the large sample of eROSITA sources will require extensive multiwavelength follow-up observations.

Despite these difficulties, the recent release of the 4XMM-DR10 catalogue, containing more than half a million X-ray sources seen with XMM–Newton, motivated us to perform a new search for INS candidates. This paper is organized as follows. In Section 2, we describe how we characterized the X-ray spectrum of the INS in a suitable way to browse the XMM–Newton source catalogue. We then describe the filtering process to exclude extended and variable sources, and those with a bright optical and/or IR counterpart (Section 3), and, finally, we analyse the new possible INSs candidates (Section 4). The results are discussed in Section 5.

2 X-RAY HARDNESS RATIOS OF 4XMM-DR10 SOURCES

The 4XMM-DR10 source catalogue was released on 2020 December 10 (Webb et al. 2020). It contains 849 991 source detections drawn from 11 647 XMM–Newton EPIC observations obtained in the first 20 yr of satellite operations. The observations cover ∼1192 deg² of the sky and contain 575 158 unique sources.

For each detection, the catalogue lists many parameters derived by the pipelines of the Science Analysis System (SAS).¹ These include the sky coordinates, the extension of the source, some flags indicating the quality of the detection and the possible variability, and the net (background subtracted) count rates in five energy bands: (1) 0.2–0.5 keV, (2) 0.5–1 keV, (3) 1–2 keV, (4) 2–4.5 keV, and (5) 4.5–12 keV.

Fig. 1 shows the HR₁ and HR₂ values of 4XMM-DR10 point-like sources (only sources with HR errors ≤0.1 are plotted). The majority of the sources are clearly grouped in two distinct regions. It has been shown in previous works (see e.g. Lin, Webb & Barret 2012) that the active galactic nuclei (AGNs, comprising BL Lac, Blazars, Quasars...) are located in the centre of the HR₁–HR₂ plane, while the lower right corner is mostly populated by stars. The blue line corresponds to the expected HR values (see below) for power-law spectra of photon index Γ = 4 and different values of absorption. AGNs, which typically have Γ < 4, lie to the left of the blue line. The spectra of stars can be described by hot plasma thermal models. The red line corresponds to emission from a collisionally ionized gas with temperature kT = 1.05 keV and different absorption values. Most stars have a lower temperature and lie below this line.

Figure 1.

HR₁–HR₂ plot of the point-like 4XMM-DR10 sources with errors on the hardness ratios smaller than 0.1 ((the error bars in the top left-hand corner illustrating this uncertainty). The sources cluster into two main groups: The AGNs, having a non-thermal spectrum, are placed to the left of the blue line (powerlaw, Γ ≲ 4); the stars, having a thermal spectrum, are placed in the region below the red line indicated by the arrow (apec, kT ≲ 1.05 keV, 0.5 solar abundance). The green line represents a blackbody (bbodyrad) of temperature kT = 0.05 keV; typical INSs and super-soft sources (SSSs), with temperatures in the range of 0.05–0.25 keV, can be found in the region above this line indicated by the green arrow. The lines are obtained by varying N_H from 0 to 10²³ cm⁻² using the tbabs model.

Open in new tabDownload slide

Compact objects form a third group of sources, comprising INSs, supernova remnants (SNRs) (that appear point-like when located in other galaxies), and binary systems hosting a white dwarf (WD), a neutron star or a stellar mass black hole. Their spectral and timing properties can be used to distinguish several classes (Fabbiano 2006): The SSSs, cataclysmic variables (CVs), novae, high-mass and low-mass X-ray binaries (HMXBs and LMXBs), and ultraluminous X-ray sources (ULXs). The group of compact objects does not have a specific location on the HR₁–HR₂ plane. However, in the lower left-hand corner of the plane (HR₁ ≲ 0 and HR₂ ≲ −0.5), there are only cold INSs and SSSs (the green line represents a blackbody with kT = 0.05 keV, and hotter blackbodies are above it).

Knowing the spectral response of the X-ray detector, it is possible to compute the expected HR values for any specific source emission model and value of the interstellar absorption. For single-component models, one can fix the parameter (e.g. the photon index or the temperature) and vary the absorption to obtain tracks on the HR₁–HR₂ plane, such as those shown by the lines in Fig. 1.

In order to derive HR values appropriate to select thermally emitting INSs, we considered a blackbody model with temperature kT ∈ [0.05, 0.25] keV and absorption column density N_H  ∈ [10²⁰, 10²²] cm⁻². We used the models implemented in xspec, and the interstellar absorption with the Tuebingen-Boulder ISM absorption cross-section, with abundances from Wilms, Allen & McCray (2000). Folding the model with the response functions of the EPIC-pn camera for different operating modes, optical filters, and off-axis angles, we computed the expected count rates in the five standard energy bands defined above. We found significant differences in the resulting HR values only as a function of the optical filters. In fact, the HRs computed for different observing modes and off-axis angles differ by less than 0.05. Therefore, in the following, we adopted, for each of the 3 filters (thin, medium, and thick) the corresponding HR values computed with the on-axis full-frame response matrices.

Fig. 2 shows the HR₁–HR₂ plane for the thin filter. It is clear, that, for the adopted model and range of parameters, only a specific region is occupied (HR₂ < HR₁), and that to each (kT, N_H) pair corresponds a (HR₁, HR₂) pair. The temperature increases from low HR₂ to high HR₂ (see the left-hand panel of Fig. 2), while the absorption increases from low HR₁ to high HR₁ (see the right-hand panel of Fig. 2).

Figure 2.

HR₁–HR₂ plane of EPIC-pn sources computed with full-frame operating mode, thin filter, and on-axis. The left-hand panel shows how HRs change varying the blackbody temperature from 0.05 to 0.25 keV; The right-hand panel shows how HRs change varying the absorption from 10²⁰ to 10²² cm⁻².

Open in new tabDownload slide

We performed polynomial fits to the curves that include the interesting region characterizing soft sources: The curves delimiting the left-hand and the right-hand sides were obtained by varying kT and keeping N_H fixed to 10²⁰ and 10^22¬cm⁻², respectively; those for the lower and the upper sides were obtained by varying N_H and keeping kT fixed to 0.05 and 0.25 keV, respectively. We fitted the left-hand and the lower sides with a function |$f_j (\mathrm{HR}_1)=\sum _i c_{ij}\mathrm{HR}_1^i$|⁠, while the right-hand and the upper sides with |$f_j (\mathrm{HR}_2)=\sum _i c_{ij}\mathrm{HR}_2^i$|⁠, where j refers to the different side and c_ij is the polynomial coefficient with i varying from 0 to 6 at maximum. All the derived coefficients are listed in Table 1.

We performed the same kind of analysis on the HR₃–HR₄ plane, and we found that the numerical HRs are restricted in the lower left corner of the plane. Considering an error of 0.05, we obtained consistent values of maximum HRs for the three filters: |$\mathrm{HR}_{3,\rm max} = -0.78(5)$| and |$\mathrm{HR}_{4,\rm max} = -1.00(5)$|⁠.

3 SELECTION OF INS CANDIDATES

Among all 849 991 detections contained in the 4XMM-DR10 catalogue, we first excluded spurious detections (sc_sum_flag < 4) and those observations in which the pn was not operating. This reduced the total sample to 688 081 detections, corresponding to 496 716 unique sources. Then, we excluded 60 914 sources that were flagged as extended (sc_ext_ml< 4) or as variable within the single observations (sc_var_flag ≠ ‘T’).

In order to obtain reliable estimates of the source spectral shape, among the remaining sources, we retained only the detections with small HR errors, namely pn_hr1_err ≤ 0.1 and pn_hr2_err ≤ 0.1. This significantly reduced the number of selected detections: 34 141, corresponding to 24 961 unique sources.

We searched for optical and/or infrared (IR) counterparts of the selected sources in the following catalogues: USNO A2.0 (Monet 1998), USNO B1.0 (Monet et al. 2003), GAIA DR2 (Gaia Collaboration et al. 2018), SDSS DR12 (Alam et al. 2015), Pan-STARRS1 (Skrutskie et al. 2016), and 2MASS (Cutri et al. 2003).

We also considered the work of Khan (2017), which catalogues point-like sources in the vicinity of nearby galaxies (<1.9 Mpc) observed by the Spitzer Space Telescope (3.6–8 and 24 μm, Werner et al. 2004). This allowed us to identify mid-IR luminous stars in the crowded and dusty discs of large star-forming galaxies, which are not detected in the near-IR and optical bands.

The cross-correlation task was performed using the online CDS X-Match Service.² The 4XMM-DR10 catalogue provides for each source the radius of the statistical error region and a systematic error, that are added in quadrature to give the total positional uncertainty. Based on this uncertainty and on the positional errors of the used catalogues, we adopted a threshold in the correlation radius corresponding to a significance of 3σ, following the prescription of Pineau et al. (2011).

As a final step of the filtering process, we wanted to select the softest sources contained in the HR region, as shown in Fig. 2, where for each observation we used the boundaries corresponding to the appropriate filter (see Table 1). We added in quadrature a systematic value of 0.05 to the statistical errors on HR to take into account the uncertainties in the predicted HR values described above. Sources were retained if their HR₁ and HR₂ were inside the boundaries and HR₃ < HR|$_{3,\rm max}=-0.78(5)$| considering their 1σ errors.

We did not use HR₄ because such soft sources are not detected or have a very low signal to noise ratio above 2 keV. In this way, we obtained 469 detections, corresponding to 140 unique sources. In Table 2, we summarize the steps of our filtering process of the 4XMM-DR10 catalogue.

4 CHARACTERIZATION OF SOFT SOURCES

All the sources matching our selection criteria were checked individually in order to verify their nature. We found that a few tens of them were actually spurious detections caused by straylight from bright sources outside the field of view, or by point-like knots of extended sources, such as bright SNRs. For other sources, we found plausible optical counterparts that were missed in the automatic cross-correlation due to a significant proper motion. We finally removed a couple of sources that lie within the soft region of the HR₁–HR₂ plane only in one detection out of many, but had average values of HR consistently outside of the region.

This further screening reduced the sample to 47 sources, six of which are not associated with any known X-ray object (see Table 3). The latter, being very soft, non-variable point-like sources, can be considered as potential INS candidates.

Many of the sources already identified are superimposed, or close to, nearby galaxies. To verify a possible association, we compared their position with the apparent dimensions of the galaxy, defined by the minor and major isophotal (m_B = 25) diameters, and their position angle, as reported in The Third Reference Catalogue of Bright Galaxies (de Vaucouleurs et al. 1991).³ Most of the sources with an associated galaxy in Table 3, are located inside the m_B = 25 ellipse. Only two sources are located at much larger distances, at ∼3 and ∼5 times the isophotal radius. The first one, 4XMM J040325.2−431721, is probably a background AGN (Ducci et al. 2014), while 4XMM J031722.7−663704 is one of the unidentified sources (see Section 4.3).

Another property that can help to classify X-ray sources is the evolution of the spectrum through different epochs. Most of the 47 selected sources were observed several times by XMM–Newton. We divided these sources into two groups: The spectrally non-variable ones, which have consistent HR values across all observations, and the spectrally variable sources, which have significantly different HR values in different epochs. This spectral variability can be visualized by tracks in the HR₁–HR₂ plane. Fig. 3 shows the HR₁ and HR₂ values of the spectrally variable (left-hand panel) and non-variable (right-hand panel) sources.

Figure 3.

The HR₁–HR₂ plane of the 47 soft sources divided into spectrally variable (left-hand panel) and non-variable (right-hand panel). They are numbered according to Table 3, and coloured according to their class: INSs (black stars), LMXBs and HMXBs (orange triangles), SNRs (blue small dots), SSSs (magenta squares), TDEs (grey diamonds), AGNs (cyan empty diamonds), and ULXs (green empty dots). The sources of unknown classification are represented by red dots.

Open in new tabDownload slide

We immediately notice the high correlation between the source class and its spectral variability: Among the variable sources, there are extra-galactic objects such as SSSs (magenta squares), ULXs (green empty dots), LMXBs (orange triangles), and one peculiar TDE, already noticed for its remarkably soft spectrum by Lin et al. (2018). We also included in the plot three sources (number 14, 15, and 20) even if they have only one XMM–Newton detection, because their variability was measured by Chandra and Swift-XRT. The first source, 4XMM J063045.4−603131, is most likely either a TDE (Mainetti, Campana & Colpi 2016) or a Galactic Nova (Oliveira et al. 2017), while the other two are LMXBs in the ω Centauri globular cluster (Cool et al. 2013) and in M31 (Barnard et al. 2014).

Strong spectral variability is also visible in the magnetar SGR J0418+5729 (Guillot et al. 2015). All the other INSs (black stars) have steady HRs: We found the CCO 1E 1207.4−5209 (Bignami et al. 2003), the high magnetic field pulsar J0726−2612 (Rigoselli et al. 2019a), two thermally emitting pulsars (McGowan et al. 2006; Ng et al. 2007), the INS candidate of Pires et al. (2015) and, as expected, six of the seven known XDINSs (the lack of RX J0420.0−5022 is due to the fact that its spectrum is so soft that its HR₂ has an error of >0.1). We also found PSR J1400−1431, a binary system composed of a millisecond pulsar (MSP) and a WD (Swiggum et al. 2017).

The other bulk of spectrally non-variable sources is made up of SNRs (blue small dots) in nearby galaxies, as M31 and M33, the AGN candidate (cyan empty diamond) previously mentioned and one HMXB in NGC 7793 (Mineo, Gilfanov & Sunyaev 2012). Finally, we showed with red dots the six INS candidates.

4.1 Analysis of the six unknown sources

To study in more detail the six sources of unknown nature and located far from nearby galaxies that we consider as potential INS candidates, we analysed other observations in the Chandra, Swift-XRT, and XMM–Newton public archives.

In order to maximize the signal-to-noise ratio for these generally faint sources, we extracted their EPIC-pn spectra using the maximum likelihood method, as described in Rigoselli et al. (2019b), Rigoselli et al. (2021). The energy bins are chosen in such a way to have a significant detection, which means at least 50 counts per bin. Spectral fitting was done with the xspec software.

For all the sources discussed below, power-law models could be rejected with high confidence either because of reduced χ² ≥ 2 or because of unphysically large photon indices (Γ ≥ 6).

4.2 4XMM J022141.5−735632

The source 4XMM J022141.5−735632 was detected only once by XMM–Newton on 2012 February 10 (obs. ID 0674110401), in an observation of the Magellanic Bridge, an H i gaseous structure connecting the Small and the Large Magellanic Clouds (Hindman, Kerr & McGee 1963).

The fit with an absorbed blackbody is in reasonably good agreement with the data (⁠|$\chi ^2_\nu =1.44$| for 16 dof, nhp = 0.11), yielding temperature kT = 0.062 ± 0.004 keV, unabsorbed 0.2–2 keV flux of |$(4.0_{-1.0}^{+1.5})\times 10^{-13}$| erg s⁻¹cm⁻², and N_H = (4 ± 1) × 10²⁰cm⁻². This N_H value is significantly smaller than the total Galactic column density in this direction, N_H,tot =1.4 × 10²¹ cm⁻² (HI4PI Collaboration et al. 2016), suggesting that this source is a nearby Galactic object. For an assumed distance of 1 kpc, the blackbody emitting radius is  |$R=6.7_{-1.7}^{+2.5}$| km and the luminosity is |$(4.8_{-1.2}^{+1.8})\times 10^{31}$| erg s⁻¹. Other fits with single component models gave worse χ² values.

A better fit can be obtained by adding to the blackbody a second spectral component, such as a power law or a hotter blackbody, but, owing to the small number of counts above 1 keV, we had to fix the photon index or the second temperature to reasonable values in order to constrain the other parameters. In the case of the addition of a power law with photon index Γ = 3, the blackbody has kT = 0.056 ± 0.004 keV, |$R=10.8_{-3.3}^{+5.3}$| km, and |$F = (5.9_{-1.8}^{+2.9})\times 10^{-13}$| erg s⁻¹cm⁻² (⁠|$\chi ^2_\nu = 0.9$| for 15 dof). In the case of a hot blackbody with kT_hot = 0.15 keV, we obtained kT = 0.049 ± 0.004 keV, |$R=21_{-7}^{+13}$| km, and |$F = (10_{-4}^{+7})\times 10^{-13}$| erg s⁻¹ cm⁻² (⁠|$\chi ^2_\nu = 0.6$| for 15 dof). In both cases, the fitted absorption column is smaller than N_H,tot and the spectral parameters are plausible for a INS at ∼1 kpc.

On the other hand, if we assume that the source is in the Magellanic Bridge, its luminosity would be of the order of 10³⁵ (d/60 kpc)² erg s⁻¹. Although this luminosity is consistent with that of an X-ray binary, the absence of an optical counterpart and lack of variability (see below) disfavour this possibility.

4XMM J022141.5−735632 was in the field of view of Swift-XRT seven times, between 2011 and 2015. The individual measurements are consistent with a constant count rate of (2.23 ± 0.55) × 10⁻³ cts s⁻¹, which assuming the best-fitting parameters of the single-blackbody fit, corresponds to a flux of F = (2.3 ± 0.6) × 10⁻¹³. The light curve shown in Fig. 4 does not give evidence for significant variability on long time-scales.

Figure 4.

Long-term light curves of the sources with multiple observations. The flux, corrected for absorption, is evaluated between 0.2 and 2 keV. Stars: values measured from a spectral fit; dots: values derived with the FTOOL pimms with the best-fitting parameters; black points: XMM–Newton EPIC-pn; red points: Chandra ACIS; blue points: Swift-XRT. Top left panel: 4XMM J022141.5−735632; top right-hand panel: 4XMM J175437.8−294149; middle left-hand panel: 4XMM J181844.3−120751; middle right-hand panel: 4XMM J220221.4+015330; bottom panel: 4XMM J031722.7−663704.

Open in new tabDownload slide

4.3 4XMM J031722.7−663704

4XMM J031722.7−663704 is the softest source of our sample, with almost no signal above 1 keV. It was observed many times by XMM–Newton, Chandra, and Swift-XRT because it is located close to the galaxy NGC 1313, which hosts two ULXs (Bachetti et al. 2013) and the ultraluminous supernova SN 1978K (Ryder et al. 1993). The angular distance of ∼9 arcmin from the centre of the galaxy is about five times larger than the isophotal radius, implying that 4XMM J031722.7−663704 is likely not associated with NGC 1313, as also noticed by Liu (2011).

Only 10⁴ of the about 30 available XMM–Newton observations that span more than 17 yr, have enough counts to perform a spectral analysis. We fitted these spectra simultaneously with an absorbed blackbody, linking the column density to a common value and obtaining |$\chi ^2_\nu =0.80$| for 58 dof (nhp = 0.86). The best-fitting N_H is (1.3 ± 0.2) × 10²¹ cm⁻², while temperatures and emitting radii were consistent within the errors, with average values of kT ≈ 0.06 keV and R ≈ 7 km. The fluxes, represented in the lowermost panel of Fig. 4 by black stars, are consistent with a constant value of ∼1.2 × 10⁻¹³erg s⁻¹ cm⁻² (⁠|$\chi ^2_\nu =1.17$| for a fit with a constant).

Imposing a common kT gave a similar |$\chi ^2_\nu =0.89$| for 67 dof [N_H = (1.3 ± 0.2) × 10²¹ cm⁻² and kT = 0.060 ± 0.003], while imposing the same parameters in all the spectra gave a worse fit (⁠|$\chi ^2_\nu =1.48$| for 76 dof). The spectral parameters are summarized in Table 4.

The remaining XMM–Newton data, for which spectral analysis was not feasible, gave detections or upper limits consistent with a constant flux (see Fig. 4). We also checked that the observations obtained with Chandra (3.6 ± 2.0 net counts on 2003 October 02, obs. ID 3551) and Swift-XRT (hundreds of observations from 2006 February 3 to 2021 April 25) do not give evidence for variability.

4.4 4XMM J175437.8−294149

The source 4XMM J175437.8−294149 was detected four times by XMM–Newton and twice by Chandra, in the course of a campaign started in 2004 to characterize the X-ray sources of the Galactic bulge in a region of low extinction, called ‘Stanek’s window’ (van den Berg et al. 2006). Only one Chandra (obs. ID 4547) and one XMM–Newton (obs. ID 0402280101) observations had enough counts for a spectral analysis.

We fitted the two spectra with an absorbed blackbody, imposing a common absorption value but letting the other parameters free to vary. We found a good fit (⁠|$\chi ^2_\nu =0.72$| for 8 dof) with N_H|$= (2.5_{-1.2}^{+1.7})\times 10^{21}$| cm⁻², temperatures kT_Chandra = 0.107 ± 0.015 keV, kT_XMM = 0.120 ± 0.025 keV, and fluxes |$F_{\rm Chandra} = (1.3_{-0.7}^{+2.7})\times 10^{-13}$| erg s⁻¹cm⁻², |$F_{\rm XMM} = (1.5_{-0.9}^{+4.9})\times 10^{-13}$| erg s⁻¹cm⁻² (unabsorbed, 0.2–2 keV).

Given that temperatures and fluxes were consistent with the same values, we fitted the two spectra linking all the parameters. This gave N_H|$= (3.2_{-1.1}^{+1.4})\times 10^{21}$| cm⁻², kT = 0.102 ± 0.014 keV, and |$F=(1.1_{-0.5}^{+1.1})\times 10^{-13}$| erg s⁻¹ cm⁻² (⁠|$\chi ^2_\nu =0.85$| for 10 dof). We used these best-fitting parameters, also reported in Table 4, to convert the count rates of the other Chandra (obs. ID 5303) and XMM–Newton detections (obs. IDs 0206590101, 0206590201, and 0801683001) into unabsorbed fluxes. The source flux remained constant from 2004 to 2018 (see Fig. 4).

The best-fitting N_H is similar to the total column density in the direction of the source, N_H,tot = 3.9 × 10²¹cm⁻² (HI4PI Collaboration et al. 2016), but its relatively large uncertainty does not allow us to estimate the source distance. For a reference distance of 1 kpc, the blackbody normalization corresponds to an emitting radius of ≲4.6 km and a bolometric luminosity of ≲ 10³²erg s⁻¹. This is consistent with a cooling INS at distances up to a few kpc.

4.5 4XMM J180528.2−273158

The source 4XMM J180528.2−273158 was detected in an observation (obs. ID 0305970101) of the Galactic Centre region (l = |$0{_{.}^{\circ}}35$|⁠, b = −|$2{_{.}^{\circ}}07$|⁠). There are no X-ray observations from Swift-XRT or Chandra of this sky region.

The source spectrum can be well fitted (⁠|$\chi ^2_\nu =0.92$| for 17 dof, nhp = 0.55) by an absorbed blackbody with kT = 0.24 ± 0.02 keV and N_H = (2.6 ± 0.5) × 10²¹ cm⁻²(total absorption in this direction is N_H,tot = 2.9 × 10²¹cm⁻², HI4PI Collaboration et al. 2016). The 0.2–2 keV flux, corrected for absorption, is (1.1 ± 0.2) × 10⁻¹³erg s⁻¹cm⁻². If we assume a distance in the range 1–10 kpc, the emitting radius would be R ≈ 0.2−2 km, smaller than the size of a neutron star, but consistent with emission from a hotspot on the star surface.

4.6 4XMM J181844.3−120751

The source 4XMM J181844.3−120751 was detected by XMM–Newton in 2002 and 2014 in three observations (obs. IDs 0008820301, 0008820601, and 0740990101) aimed at characterizing two colliding wind binaries in the open cluster NGC 6604, which is at a distance of 1.7 kpc (De Becker et al. 2004, 2005; De Becker 2015). This region was never observed by Chandra or Swift-XRT.

4XMM J181844.3−120751, located at 12 arcmin from the centre of the cluster, was also noticed by De Becker et al. (2005) and Lin et al. (2012), who suggested that it is a compact object, based on its soft spectrum and the absence of an optical counterpart.

None of the three spectra gave an acceptable fit with the blackbody model (⁠|$\chi ^2_\nu =2.9$| for 19 dof, nhp = 3 × 10⁻⁵). Also, fits with other simple single-component models (power law, bremsstrahlung, thermal disk), as well as with the combination of a blackbody plus power-law, were rejected (⁠|$\chi ^2_\nu \gt 3$|⁠).

A good fit (⁠|$\chi ^2_\nu =1.2$| for 19 dof, nhp = 0.22) was instead obtained with a model of collisionally ionized diffuse gas with abundances fixed at solar values (apec in xspec). The best-fit N_H is |$(5.5_{-0.6}^{+0.7})\times 10^{21}$| cm⁻², similar to N_H,tot =7.1 × 10²¹ cm⁻² (HI4PI Collaboration et al. 2016) and to the values obtained for the two binaries in NGC 6604 (De Becker et al. 2005; De Becker 2015).

The three temperatures are consistent within errors (kT₁ = 0.29 ± 0.03 keV, kT₂ = 0.33 ± 0.04 keV, and kT₃ = 0.27 ± 0.03 keV), as well as the unabsorbed fluxes (see Fig. 4). Imposing the same parameters in all the spectra gave a still acceptable fit (⁠|$\chi ^2_\nu =1.25$| for 23 dof), with kT = 0.30 ± 0.03 keV and |$F=(1.2_{-0.4}^{+0.8})\times 10^{-12}$| erg s⁻¹cm⁻².

4.7 4XMM J220221.4+015330

4XMM J220221.4+015330, the brightest source of our sample, was detected three times in Swift-XRT (obs. ID 00036084001), Chandra (obs. ID 10142), and XMM–Newton (obs. ID 0655346833) short observations. The three corresponding spectra have a good signal-to-noise ratio and clearly indicate spectral variability. They could not be fitted by either a blackbody or a power law (⁠|$\chi ^2_\nu \gt 2$|⁠). Acceptable fits were obtained by adding a second component, but only when the cooler component has a normalization free to vary among the epochs.

The power-law plus blackbody best-fitting parameters are N_H = (5 ± 2) × 10²⁰ cm⁻², Γ = 3.0 ± 0.3, N = (1.5 ± 0.5) × 10⁻⁴photons cm⁻²s⁻¹keV⁻¹ at 1 keV and kT = 0.188 ± 0.009 keV (⁠|$\chi ^2_\nu =1.23$| for 40 dof, nhp = 0.15). The blackbody flux shows significant variations across the three observations: F_XRT = (5.6 ± 1.5) × 10⁻¹³ erg s⁻¹ cm⁻², F_Chandra = (4.0 ± 0.9) × 10⁻¹³erg s⁻¹ cm⁻², and F_XMM = (9.1 ± 0.8) × 10⁻¹³ erg s⁻¹ cm⁻².

The two-blackbody model gives a better |$\chi ^2_\nu =1.10$| (nhp = 0.31), a small absorption column (N_H < 2.7 × 10²⁰cm⁻²), |$kT_{\rm hot}=0.33_{-0.03}^{+0.04}$| keV, and kT = 0.16 ± 0.01 keV. The hotter blackbody has a radius of |$0.16_{-0.05}^{+0.06}$| km for a distance of 1 kpc. The long-term light curve of the colder blackbody is shown in Fig. 4. The fluxes of the cooler blackbody vary similarly to the previous case: F_XRT = (6.6 ± 1.5) × 10⁻¹³erg s⁻¹cm⁻², |$F_{\rm Chandra}=(4.6_{-1.0}^{+0.9})\times 10^{-13}$| erg s⁻¹ cm⁻², and |$F_{\rm XMM}=(9.8_{-1.1}^{+1.2})\times 10^{-13}$| erg s⁻¹ cm⁻²).

We notice that a good fit (⁠|$\chi ^2_\nu =1.17$| for 42 dof, nhp = 0.21) can also be achieved with an absorbed bremsstrahlung, yielding N_H = (7 ± 1) × 10²⁰ cm⁻² and kT = 0.53 ± 0.03 keV. The fluxes are higher than the previous ones, but a significant variability is still present: F_XRT = (1.7 ± 0.2) × 10⁻¹² erg s⁻¹ cm⁻², F_Chandra = (1.2 ± 0.1) × 10⁻¹²erg s⁻¹cm⁻², and F_XMM = (2.3 ± 0.2) × 10⁻¹²erg s⁻¹ cm⁻².

The source has a possible optical counterpart with magnitudes R = 25.1 and K_s = 22.5 detected by VIPERS-MLS (Moutard et al. 2016). From its high X-ray flux ∼10⁻¹² erg s⁻¹cm⁻², we derived log (F_X/F_O) = 3.4 and log (F_X/F_IR) = 4.0, which are typical of a compact object. However, the observed variability excludes 4XMM J220221.4+015330 from the list of the INS candidates.

5 DISCUSSION AND CONCLUSIONS

We searched for new INS candidates in the 4XMM-DR10 source catalogue, using selection criteria based on their spectral shape, as inferred from X-ray HRs in the softest energy bands (0.2–0.5, 0.5–1, and 1–2 keV), and on cross-correlations with catalogues of possible optical/IR counterparts. In particular, we identified a region in the HR₁–HR₂ plane, assuming a blackbody spectrum absorbed by the interstellar medium folded through the response of the EPIC-pn detector, and used it to select the softest sources.

From more than half a million X-ray sources contained in the 4XMM-DR10 catalogue, we finally obtained a sample of 47 point-like sources. This sample includes about twenty SNRs or X-ray binaries located in nearby galaxies, a few AGNs and, as expected, several already-known INSs. The remaining six sources have an unknown nature and are located far from nearby galaxies, therefore we considered them as potential INS candidates. A more detailed spectral and timing analysis, using also Chandra and Swift-XRT data, showed that two of them are unlikely to be INSs. 4XMM J181844.3−120751 has a spectrum inconsistent with blackbody emission, and is instead well fit with a thermal plasma model with kT ≈ 0.3 keV. This source is close to the open cluster NGC 6604, but the lack of a bright optical counterpart disfavours a stellar origin. 4XMM J220221.4+015330 showed clear spectral and flux variability. Its faint optical counterpart implies log (F_X/F_O) = 3.4, consistent with a peculiarly soft AGN or an X-ray binary.

The remaining four sources (4XMM J022141.5−735632, 4XMM J031722.7−663704, 4XMM J175437.8−294149, and 4XMM J180528.2−273158) have soft thermal spectra and show no evidence for long-term variability. Their temperatures and emission radii, as inferred from blackbody fits, are consistent with emission from hot spots or from (a large fraction of) the whole surface of INSs. The spectrum of the source with the highest signal to noise ratio, 4XMM J022141.5−735632, was better fitted with two components models: Either two blackbodies (kT_hot = 0.15 keV and kT ≈ 0.049 keV), or a power law plus blackbody (Γ = 3 and kT ≈ 0.056 keV), as typically observed in INSs.

In Table 4, we also report the expected R magnitude of these candidates, obtained extrapolating their best-fitting X-ray blackbody. These magnitudes are rather faint, but we note that these can be considered as lower limits, because the optical counterparts identified for INSs are a factor ∼5–10 times brighter than the extrapolation of the Planckian spectra (Mignani 2009). Whether this is due to a deviation from blackbody emission and/or a non-uniform temperature distribution, or a contamination from a nearby diffuse source, it is still an open issue (Wang et al. 2017, 2018). Deep multiwavelength observations are needed to identify the counterparts of these four sources, and possibly confirm the suggested neutron star nature, e.g. through the measurement of proper motions.

The small number of candidate INSs we found in the 4XMM-DR10 catalogue is not surprising, considering the rarity of thermally emitting neutron stars sufficiently close to be detected in soft X-rays (Posselt et al. 2008). The LogN–LogS derived by these authors with a detailed evolutionary synthesis model, and taking into account the nonuniform distribution of interstellar absorption, predicts the presence of the order of ∼100 cooling INSs in the whole sky above 0.001 ROSAT counts  s⁻¹ (corresponding to a 0.2–2 keV flux of a few 10⁻¹⁵erg s⁻¹ cm⁻²). The XMM–Newton observations used in the 4XMM-DR10 catalogue have a non uniform sensitivity, down to F_0.2−2 ∼ 10⁻¹⁶ erg s⁻¹cm⁻² for the faintest detected sources. However, our requirement of having HRs with errors <0.1 limited our search to sources with fluxes above ∼10⁻¹⁵ erg s⁻¹cm⁻². Considering that the sky coverage of the catalogue is only ∼3 per cent, the small number of candidates found in our work is consistent with the expectations. It is possible that other INSs are among the weakest 4XMM-DR10 sources for which an accurate spectral characterization is currently impossible.

ACKNOWLEDGEMENTS

The scientific results reported in this article are based on data obtained from XMM–Newton, Swift, and Chandra. The data analysis has benefited from the software provided by the NASA’s High Energy Astrophysics Science Archive Research Center (HEASARC). We extensively made use of the Simbad, VizieR and X-Match databases, operated at CDS, Strasbourg, France. We acknowledge the whole XMM–Newton SSC Consortium for the 4XMM-DR10 catalogue production. We thank F. Pintore for help with the Chandra data and the referee F. Walter for useful comments. We acknowledge financial support from the Italian Ministry for University and Research through grant 2017LJ39LM ‘UNIAM’ and the INAF ‘Main-streams’ funding grant (DP n.43/18).

DATA AVAILABILITY

All the data used in this article are available in public archives.

Footnotes

REFERENCES