© A. M. Pires et al. 2022

1. Introduction

The vast majority of the neutron stars in our Galaxy are observed at radio wavelengths. Yet, it is arguably in X-rays that isolated neutron stars (INSs) reveal their diversity in all their complexity (see e.g. Kaspi 2010, for an overview). In particular, the group of radio-quiet thermally emitting INSs dubbed the “magnificent seven” (M7, for simplicity), originally identified in the ROSAT All-Sky Survey data as soft, bright X-ray sources (observed fluxes typically within ∼10⁻¹² − 10⁻¹¹ erg s⁻¹ cm⁻²) with no obvious optical counterparts, share a rather well-defined set of properties that have never been encountered together in previously known families of INSs (see Haberl 2007; Turolla et al. 2009, for reviews). The M7 are locally as numerous as young radio and γ-ray pulsars and may belong to a formerly neglected component of the overall Galactic INS population. The discovery of similar sources beyond our local volume is therefore key to understanding their properties as a group and relations to other families of Galactic INSs (e.g. Popov et al. 2010; Viganò et al. 2013).

At fainter X-ray fluxes (f_X ≲ 1.5 × 10⁻¹³ erg s⁻¹ cm⁻²), source confusion and contamination from other classes of X-ray emitters hamper the identification of new members due to the large positional and spectral uncertainties. In preparation for the full sensitivity of the eROSITA All-Sky Survey (eRASS; Predehl et al. 2021), the use of serendipitous data from the XMM-Newton Observatory (Webb et al. 2020) provides an excellent opportunity to test search algorithms and discover new INSs beyond the solar vicinity.

Building on our experience cross-correlating previous releases of the XMM-Newton catalogue of X-ray sources (Pires et al. 2009; Motch et al. 2017), we searched the source content of 4XMM-DR9 for new INS candidates. We used as main criteria the absence of any catalogued optical/UV/IR counterpart and a soft spectrum down to a limiting flux of 10⁻¹⁴ erg s⁻¹ cm⁻² in the 0.5 − 1 keV energy band (see Appendix A, for details). Among the selected sources, we retrieved known thermally emitting INSs: the M7, quiescent neutron stars in globular clusters, and our “own” cooling INS identified in the Carina Nebula, 2XMM J104608.7−594306 (see Pires et al. 2015, and references therein). The brightest among the unknown sources were then selected for dedicated observations with XMM-Newton (fulfil programs 088419, 090126; Pires et al., in prep.).

We report the first results of this follow-up effort on 4XMM J022141.5−735632, an X-ray source located in the direction of the Magellanic Bridge and independently identified by similar searches in 4XMM-DR10 and eRASS data (Rigoselli et al. 2022b; Kurpas et al., priv. comm.). The overall properties of 4XMM J022141.5−735632, as obtained from the 4XMM-DR9 catalogue, are listed in Table 1.

In the following, we present the XMM-Newton and eROSITA data set that forms the backbone of the analysis; we show that the X-ray emission of the target is predominantly thermal and constant over many years; the absence of optical counterparts in Legacy Survey¹ (Dey et al. 2019) Data Release 10 imaging safely excludes a more ordinary X-ray emitter than an INS. We finalise with a discussion of our state of knowledge and interpretation on the nature of the source, with prospects for future investigations.

2. Observations and data reduction

2.1. XMM-Newton

The XMM-Newton Observatory (Jansen et al. 2001) targeted the X-ray source 4XMM J022141.5−735632 (J0221, for short) for circa 50 ks on July 9–10 2021. We included in the analysis the only other XMM-Newton observation that serendipitously detected the target in 2012 (0674110401). We performed standard data reduction with SAS 20 (xmmsas_20211130_0941) using the most up-to-date calibration files and following the analysis guidelines of each EPIC instrument. We filtered the event lists to exclude “bad” CCD pixels and columns and retain the photon patterns with the highest quality energy calibration. The source centroid and optimal extraction region in each EPIC camera, with typical sizes of 30″ − 40″, were defined with the SAS task eregionanalyse in the 0.2 − 2 keV energy band; the X-ray emission of the source is compatible with the background level at energies above 2 keV. Background circular regions of size 100″ were defined away from the source, on the same CCD of the target whenever possible. We show in Table 2 the instrumental setup, net exposure, and percentage of good-time intervals (GTIs) of the XMM-Newton observations included in the analysis.

The parameters of the target (the detection likelihood, counts, rates, and hardness ratios, as extracted from a maximum likelihood PSF fitting on the EPIC images of each XMM-Newton epoch) are listed in the two first columns of Table 3. The parameters were determined with the SAS task edetect_chain on images created for each camera, observation, and standard catalogue energy bands (see caption of Table 1, for a definition). Only the combined EPIC results are shown.

We used the SAS task eposcorr to refine the astrometry by cross-correlating the list of EPIC X-ray source positions with those of optical Guide Star Catalogue (Lasker et al. 2021) objects lying within 15′ around the nominal pointing coordinates. The small positional offsets in right ascension and declination listed in Table 3 were consistently detected with the Gaia EDR3 and AllWISE catalogues (Gaia Collaboration 2020; Cutri et al. 2021). The agreement between epochs is good overall (< 3σ).

Stacked source detection on both XMM-Newton observations was performed using the XMM-Newton SAS task edetect_stack (Traulsen et al. 2019) and the 4XMM-DR12s catalogue pipeline. First, the results of the eposcorr astrometric correction of the individual observations were applied to the attitude information in order to shift the events to rectified positions, as described by Traulsen et al. (2020). We then extracted images, background maps, exposure maps, and detection masks for both observations with reference coordinates centred at the position of J0221. Deviating from the default parameters, we used an image binning of 2″  ×  2″ pixels and a PSF fit radius of 2′ about the sources, in order to optimise the sensitivity to faint sources in these fields, which are dominated by point-like sources and do not suffer from source confusion, and in order to achieve optimum positional accuracy of the target. Maximum-likelihood PSF source detection was run simultaneously in the five XMM-Newton energy bands. In the whole field, 199 sources are detected with a minimum detection likelihood of six in the stack or in at least one contributing observation. Figure 1 shows a false-colour mosaic image of the field, colour coding energies between 0.2 and 1 keV (red), 1 keV and 2 keV (green), and 2 keV and 12 keV (blue). For displaying purposes, it has been smoothed by a top-hat function with a two-pixel radius. The X-ray soft J0221 is prominently visible in the red band.

Fig. 1. Smoothed false-colour (top) and source detection (bottom) images of the two XMM-Newton observations covering the position of J0221. The red band covers the energies from 0.2−1.0 keV, green from 1.0−2.0 keV, and blue from 2.0−12.0 keV. The position of the target and of other sources detected in the field-of-view are marked by white and light blue circles, respectively.

2.2. eROSITA

The All-Sky Survey of eROSITA, the first to be performed at X-ray energies since the ROSAT era (Voges et al. 1999), is expected to surpass its predecessor’s sensitivity by a factor of 25 (0.2 − 2.3 keV; Predehl et al. 2021). Since the beginning of the survey in December 2019, eROSITA has successfully completed four of the planned eight all-sky scans. The eRASS catalogues are created by and made available to the members of the German eROSITA Consortium. We searched the individual (eRASS1–4) and stacked (eRASS:4) X-ray source catalogues for unidentified INSs above a limiting flux of 10⁻¹³ erg s⁻¹ cm⁻² in the 0.2 − 2 keV energy band (Kurpas et al., in prep.). Having survived our probabilistic catalogue cross-matching and screening procedure (only about 0.05% of over 10⁵ X-ray sources satisfy all the selection criteria for identification with an INS candidate), J0221 is one of a handful of intrinsically soft point sources with no obvious counterparts that have been lined up for follow-up investigations in the optical and in X-rays.

The eROSITA instrument consists of seven telescope modules (TM 1–7) operating in the 0.2 − 10 keV energy band. The detectors of TM 5 and 7 (“TM9”) are more sensitive to soft X-rays, as they lack the aluminium on-chip optical light filter the other five cameras (“TM8”) carry. The TM9 detectors are known to suffer from time-variable light leaks that affect their performance and calibration at the softest energies (see Predehl et al. 2021, for details).

For the analysis, we retrieved the event lists corresponding to the target’s sky tile in each eRASS1–4. We list in Table 2 an overview of these observations. The position of the X-ray source was repeatedly covered in April and October 2020–2021, with exposures of ∼0.9 − 1 ks. The data sets were processed with pipeline version c020 and analysed in the 0.2 − 10 keV energy band with the eSASS software system, applying up-to-date calibration files (eSASSusers_211214; Brunner et al. 2022). For all observations, the event files of each individual TM were filtered for periods of high background activity with the eSASS task flaregti in the 2.2 − 10 keV energy band. Unless otherwise noted, all valid photon patterns and active telescope modules (“TM0”) were considered for optimal sensitivity. The total sum of GTIs accumulated over the four epochs is 2.3 ks (∼1 ks, corrected for the vignetting). Altogether, 101 ± 12 photons are collected from the X-ray source in the 0.2 − 1 keV energy band, considering all seven active TMs.

Similarly to the XMM-Newton observations, we analysed the X-ray source content of the observations and the parameters of the target using maximum likelihood PSF fitting (Brunner et al. 2022). Based on the results of source detection and PSF fitting, we used the “auto” option of the eSASS task srctool to create optimised source and background extraction regions over the cumulative eRASS:4 events (Fig. 2). The positions of the 67 X-ray sources detected in the field-of-view in the 0.2 − 1 keV energy band were excluded from the background region. Finally, we verified the statistics of the eRASS light curves of the target for general trend variability. The light curves were corrected for bad pixels, vignetting, exposure, and background counts with the eSASS task srctool. The 0.2 − 1 keV count rate of the target is consistent with a constant value.

Fig. 2. eROSITA extraction regions (eRASS:4, TM0, 0.2 − 1 keV). Top panel: background annulus (magenta) with inner and outer radii of sizes 2.5′ and 15′, respectively, excluding “contaminants” (cyan circles with a red bar across). Bottom panel: source circular region (yellow) of size 82″.

3. Results

3.1. Timing analysis

We searched for periodic signals that could be associated with the rotation period of the X-ray source (hundreds of milliseconds to tens of seconds); alternatively, to identify flux modulations of the order of tens of seconds to hours as observed in novae and supersoft X-ray sources (Orio et al. 2022). We considered only the two ∼30 ks XMM-Newton observations of J0221 (Table 2), the individual eRASS exposures being too shallow for a meaningful timing analysis. For maximum sensitivity, we included all valid patterns and events of the three EPIC cameras unfiltered for GTIs in the 0.2 − 2 keV energy band. The times-of-arrival of the photons, barycentred to the rest frame of the solar system using the SAS task barycen and the source coordinates (Table 3), were Fourier-analysed in the frequency domain to search for the presence of periodic signals (Buccheri et al. 1983). We analysed the EPIC cameras together in the Δν = (0.002 − 1.9)×10⁻¹ Hz frequency range; the higher time resolution of the pn camera with respect to MOS allows one to search the pn time series in a broader frequency range, up to ∼6.8 Hz, albeit with somewhat lower sensitivity. We adopted in both (EPIC/pn) searches frequency steps of 8 − 10 μHz, which correspond to an oversampling factor of ∼3. The number of statistically independent trials are ∼(2 − 3)×10⁵ and ∼(6 − 8)×10³, respectively, for pn and EPIC searches in the 2012 and 2021 observations. We found no statistically significant (> 4σ) periodicity in neither epoch. The most constraining 3σ upper limits on the source pulsed fraction, $p_{\mathrm{f}}^{\mathrm{pn}}<13$% (P ∼ 0.15 − 5000 s) and $p_{\mathrm{f}}^{\mathrm{EPIC}}<10$% (P ∼ 5.2 − 5000 s), are derived from the 2021 observation.

3.2. Spectral analysis and long-term variability

The analysed data set comprises six epochs, ten spectra², and over 4000 counts (0.2 − 2 keV); background noise amounts to up to 15% and 30%, respectively, of the total EPIC and eROSITA events. The spectral analysis was restricted to GTI-filtered photons, single- and double-pattern events for pn, and all valid CCD patterns for the MOS and TM cameras. We regrouped the energy channels to avoid low (< 5) counts per spectral bin and kept oversampling of the instrumental resolution capped to a maximum factor of 3. To fit the spectra we used XSPEC 12.10.1f and the Cash C statistic (Arnaud 1996; Cash 1979); unless otherwise noted the fit parameters were allowed to vary freely within reasonable ranges. The spectra were fitted simultaneously, allowing for a renormalisation factor to account for cross-calibration uncertainties between epochs and instruments. We verified that the inclusion of photons from eROSITA detectors 5 and 7 (those affected by optical leaks; see Sect. 2.2) do not bias the parameter estimation of the results reported here. Finally, we adopted the photoionization cross-sections of Verner et al. (1996), the photoelectric absorption model tbabs, and elemental abundances of Wilms et al. (2000) to account for the interstellar material in the line-of-sight; while testing for variable elemental abundances, we adopted the absorption model tbvarabs and checked the improvement in the fit statistic element-wise by means of a F-test. Complimentary information can be found in Appendices B–C.

The main results of the spectral analysis are summarised in Table 4. In model (1) we assumed a single-temperature absorbed blackbody and constancy between all observations. The best solution consists of a soft blackbody with kT = 60.2 ± 1.7 eV, absorbed by a column density of $N_{\text{H}}=7.2_{-0.8}^{+0.9}\times {10}^{20}$ cm⁻²; the observed flux is f_X = 1.18(3)×10⁻¹³ erg s⁻¹ cm⁻² (0.2 − 2 keV) and the size of the emission region, at a distance of 1 kpc, is $8.7_{-1.2}^{+1.4}$ km. The inclusion of an additional component – either a second blackbody with kT = 140 − 260 eV, or a power-law tail with photon index $\mathrm{\Gamma }=4.2_{-1.3}^{+1.6}$ – is statistically significant (F-test probability of 1.8 × 10⁻⁵; see models 2 and 3 in Table 4). In both cases, there are no significant changes (> 2σ) to the dominant blackbody component, nor to N_H. In Fig. 3 we show the spectral data folded with model (2) and residuals.

Fig. 3. Best-fit double-blackbody model and residuals (fit 2 in Table 4). The two EPIC pn and four MOS data sets are colour-coded in dark grey and blue, respectively. The eROSITA spectrum (merged here for display purposes) is shown in pink.

In model (4), we allowed the parameters of model (1), except for N_H, to assume independent values in each epoch. We list in Table 4 the median and median absolute deviation of the parameters; see Table 5 and Fig. 4, for the individual values. We found no evidence of long-term variability: the deviations from a constant value are well within the expected cross-calibration uncertainties between the EPIC instruments (2% in kT and up to 3% in flux; Read et al. 2014); the source parameters agree within the 90% confidence level and the fit has a non-acceptable F-test probability with respect to model (1). The much wider spread and trends observed in eROSITA data (left panel of Fig. 4), discussed in Appendix B, may be asserted to low count statistics and remaining calibration issues at the softest end of the observatory passband.

Fig. 4. Long-term spectral and flux variation of the X-ray source. Left panel: blackbody temperature and emission radius as a function of MJD (“variable” model 4; see Tables 4–5). Right panel: long-term evolution of the 0.2 − 2 keV X-ray flux of the target, including upper limits and previous detection by other X-ray missions. The time interval extends back to the ROSAT All-Sky Survey (RASS) and pointed (PSPC) era and include data points from Swift XRT and XMM-Newton slew observations (see the text, for details). In all plots the purple horizontal shaded areas show the 1σ median absolute deviation of the parameters.

Interestingly, the analysis of the long-term evolution of the 0.2 − 2 keV X-ray flux of J0221, including upper limits and previous detection by other X-ray missions, suggests that the properties of the X-ray source are stable over 30 years (Fig. 4; right). In addition to the EPIC and eRASS observations analysed here, we employed the web-based tool High-Energy Lightcurve Generator (HILIGT³; Saxton et al. 2022) to retrieve upper limits and flux values for the source. These were determined through aperture photometry at J0221’s sky position and the Bayesian approach of Kraft et al. (1991) on archival observations from Swift XRT, ROSAT, and the XMM-Newton Slew Survey (see Ruiz et al. 2022, and references therein). We assumed an absorbed blackbody spectral model consistent with the source’s spectral shape and selected a 3σ confidence interval to estimate upper limits (single-sided 99.7% probability). The measurements are overall consistent with a constant flux.

Neutron star atmosphere models, with B = 10¹² − 10¹³ G, N_H = 1.1 × 10²¹ cm⁻² and T_eff ∼ (2.7 − 3)×10⁵ K, describe the spectrum of the source as well as the multi-component models, and can be considered an improved fit with respect to the non-magnetic case (models 5–7 in Table 4). Assuming emission from the entire stellar surface, canonical neutron star mass and radius, and a fully ionized hydrogen atmosphere (nsa in XSPEC; Pavlov et al. 1995), the model normalisation requires in all cases a nearby neutron star within 90 pc (3σ). These may be rather unrealistic assumptions for the atmosphere of non-accreting neutron stars, given the large increase in ionization potential expected under such conditions of temperature and magnetic field (see Potekhin 2014, for a comprehensive review). Other spectral models did not provide compelling results. In particular, the temperature of the best-fit partially ionized neutron star atmosphere (nsmaxg in XSPEC; Ho et al. 2008) is outside the computed model grid, log(T_eff) < 5.5, while the photon index of the best-fit power-law model is unreasonably steep (Γ ∼ 8). We note that our results are in perfect agreement with the findings of Rigoselli et al. (2022b).

3.3. Optical and near-infrared limits

We see no evidence in the optical and near-infrared for the source in the Legacy Survey Data Release 10 (DR10). The survey follows the same data reduction of Data Release 9 (Dey et al. 2019), but covers almost entirely the German part of the eROSITA sky⁴ at the depth of the Dark Energy Survey in g, r, i, z (Abbott et al. 2018). The depth and the addition of the i band increases the number of detected sources by 30% (Legacy Survey Collaboration 2022, private communication). In the infrared the surveys makes use of 8 years of NEOWISE observations (Mainzer et al. 2014). To search for the maximum limit of the available DECam archival imaging at this sky location, additional exposures in u and i were added and a depth weighted, multi-filter stack of the 32 exposures with exposure times greater than 60 s was produced by the DECam Community Pipeline (Valdes et al. 2014). Figure 5 is the resultant deep stack with the position of the INS candidate indicated. Clearly, there is no sign of a source in this blank region. The depth limit was obtained from the local measurement of the background noise (standard deviation of 0.52 data numbers), point-spread full-width at half-maximum (4.23 pixels corresponding to 1.06″), and zero-point calibration to Gaia-EDR G-band (28.68). Using the formula for the depth in an optimal Gaussian aperture with the measured values yields a G-band depth limit of 26.45 (25.89) at 3σ (5σ). These allow us to derive a 3σ lower limit on the target’s X-ray-to-optical flux ratio, log(F_X/F_V)≳3.7. We assumed a flat spectrum to translate the G-band flux to that in the V-band and adopted the best-fit F_X and N_H from model (2) in Table 4; optical de-reddening and the total extinction specific to the Gaia-EDR G-band were computed with the usual relations of Predehl & Schmitt (1995) and Cardelli et al. (1989).

Fig. 5. Photometric-weighted multi-filter stack of Legacy Survey DR10 g, r, z and DECam u, i images centred on the sky position of J0221. The blank (mG > 26.45) 3σ confidence level error circle of the X-ray source is displayed in the inset in an inverted black and white colour map.

4. Discussion and outlook

We report the first results of a campaign to follow-up INS candidates from 4XMM-DR9. The X-ray source 4XMM J022141.5−735632 was put forward by Rigoselli et al. (2022b) on the same premise of a soft spectrum and lack of obvious counterparts; it is also a “target of interest” on dedicated searches in eROSITA All-Sky Survey data. The joint analysis of the XMM-Newton and eROSITA observations performed between 2012 and 2021 confirms the source’s essentially thermal energy distribution, with kT ∼ 60 eV, N_H ∼ 7 × 10²⁰ cm⁻², and f_X ∼ 1.2 × 10⁻¹³ erg s⁻¹ cm⁻² (0.2 − 2 keV). The optical limit derived from the deep stacked Legacy Survey DR10 and additional DECam images, m_G > 26.45 (3σ) in the Gaia-EDR G-band, excludes a cataclysmic variable or hot white dwarf in the foreground of the Magellanic Bridge.

We find no evidence for variability in either flux or spectral state within the nearly ten-year interval covered by the analysis. Previous detection at a similar flux level by ROSAT PSPC and the Swift X-ray Telescope suggest that the emission is fairly stable over decade-long time scales. Additional monitoring with higher quality data will be necessary to exclude the few-percent level of spectral variation as reported for the M7 INSs RX J0720.4−3125 and RX J1605.3+3249 (de Vries et al. 2004; Pires et al. 2022).

In addition to the main thermal component, excess emission above 0.7 keV may be accommodated with either a second (hot) blackbody with $k{T}_{\mathrm{hot}}={190}_{-50}^{+70}$ eV or a power-law tail; for the latter, the photon index $\mathrm{\Gamma }=4.2_{-1.3}^{+1.6}$ is considerably steeper than what is typically observed in the spectra of middle-aged spin-powered pulsars dominated by thermal emission (e.g. Schwope et al. 2022; Rigoselli et al. 2022a). Magnetised and fully ionized neutron star atmosphere models, with B = 10¹² − 10¹³ G and T_eff = (2.6 − 3)×10⁵ K, are also in reasonable agreement with the data. All best-fit nsa models in Table 4 predict a rather close-by neutron star within 90 pc (3σ). We caution, however, against the applicability of such simple model atmospheres for sources as soft as J0221 when not in the weak-field regime (Potekhin 2014).

In Fig. 6 we show the range of luminosity and emission radius of J0221 in the context of the sample of 55 cooling INSs from Potekhin et al. (2020, we refer to their review paper for the terminology, references, and a description of the included data sets and spectral models). The distance interval of 200 pc to 2.9 kpc, indicated by the arrow in the diagram, assumes that the X-ray source has a similar emission radius as the M7. Considering a fiducial timescale of 1.2 Myr to cool down the surface to its present-day temperature of ∼(3 − 7)×10⁵ K (see e.g. Yakovlev & Pethick 2004, their Fig. 1), the location of J0221 ∼100 pc to 1.9 kpc below the Galactic disk (b = −41.7°; Table 1) implies a projected speed within 90 km s⁻¹ and 1600 km s⁻¹, roughly in agreement with the range observed for radio pulsars (Hobbs et al. 2005).

Fig. 6. Thermal luminosity vs. emission radius for cooling neutron stars. The arrow shows the distance-luminosity range of the target assuming the same emission radii observed for the M7 INSs (see the text, for details).

From all our tested spectral models we infer a hydrogen column density a factor of ∼1.5 − 3 higher than the Galactic value from Dickey & Lockman (1990, N_H = 3.8 × 10²⁰ cm⁻²), but within that derived from the HI4PI map in this line-of-sight (N_H = 1.4 × 10²¹ cm⁻²; HI4PI Collaboration 2016). The latter integrates the HI emission over both the Galactic and the Magellanic Cloud systems (see also Sasaki et al. 2022). If the absorption is purely interstellar, this would place J0221 outside the Milky Way into the Magellanic Bridge, at a distance of the order of that of the Small Magellanic Cloud, 60 kpc. The inferred size of the emission area for such a large distance is incompatible with the neutron star surface. On the other hand, part of the absorption could be caused by matter local to the X-ray source (e.g. the case of RX J0806.4−4123; Posselt et al. 2018). Interestingly, residuals around 540 − 580 eV may indicate an enhancement of oxygen local to or in the direction of the INS candidate (Appendix C). An updated value of the HI emission in this line-of-sight excluding the contribution of the Magellanic system is necessary for a reliable estimate of the amount of absorption intrinsic to the X-ray source.

X-ray emitting INSs at faint fluxes will consist of the products of more distant massive star associations and will be younger and hotter than the M7 on average (see Pires et al. 2017, for details based on a population synthesis model and a discussion of multiwavelength follow-up strategies of new INS candidates). The most compelling candidates will be located in the Galactic plane or close to open stellar clusters and supernova remnants. Other interesting Galactic and extragalactic X-ray sources (high-redshift quasars, ultraluminous and supersoft X-ray sources, millisecond and rotation-powered pulsars) may also be selected by our searches (Pires et al., in prep.; see also Khokhriakova et al. 2022; Vahdat et al. 2022). In particular, the sky location and softness of J0221 could indicate a supersoft nature (Kahabka & van den Heuvel 1997, for a review). However, at the distance of the SMC, the luminosity of ∼2 × 10³⁵ erg s⁻¹ is well below the range for stable hydrogen burning typical of white dwarves in close binary systems (Wolf et al. 2013). The size of the inferred emission region, ∼290 − 410 km, is likewise inconsistent with the white dwarf radius. Given the absence of significant variability and optical counterparts, we regard the surpersoft source interpretation as unlikely.

Overall, the identification of 4XMM J022141.5−735632 with an INS appears robust. A systematic search for pulsations from dedicated searches in X-rays and in the radio regime is crucial to establish the true nature of the X-ray source and shed light on its relations to other Galactic neutron stars.

Acknowledgments

We thank the anonymous referee for useful comments and suggestions which helped to improve the paper. This work was supported by the project XMM2ATHENA, which has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement n^o101004168. IT gratefully acknowledges the support of SSC work at AIP by Deutsches Zentrum für Luft- und Raumfahrt (DLR) through grant 50 OX 1901. DT acknowledges support by DLR grants FKZ 50 OR 2203. This research has made use of data obtained from the 4XMM XMM-Newton serendipitous source catalogue compiled by the 10 institutes of the XMM-Newton Survey Science Centre selected by ESA. This work is based on data from eROSITA the soft X-ray instrument aboard SRG, a joint Russian-German science mission supported by the Russian Space Agency (Roskosmos), in the interests of the Russian Academy of Sciences represented by its Space Research Institute (IKI), and the Deutsches Zentrum für Luft- und Raumfahrt (DLR). The SRG spacecraft was built by Lavochkin Association (NPOL) and its subcontractors, and is operated by NPOL with support from the Max Planck Institute for Extraterrestrial Physics (MPE). The development and construction of the eROSITA X-ray instrument was led by MPE, with contributions from the Dr. Karl Remeis Observatory Bamberg & ECAP (FAU Erlangen-Nuernberg), the University of Hamburg Observatory, the Leibniz Institute for Astrophysics Potsdam (AIP), and the Institute for Astronomy and Astrophysics of the University of Tübingen, with the support of DLR and the Max Planck Society. The Argelander Institute for Astronomy of the University of Bonn and the Ludwig Maximilians Universität Munich also participated in the science preparation for eROSITA. The eROSITA data shown here were processed using the eSASS/NRTA software system developed by the German eROSITA consortium. The Legacy Surveys consist of three individual and complementary projects: the Dark Energy Camera Legacy Survey (DECaLS; Proposal ID #2014B-0404; PIs: David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO Prop. ID #2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band Legacy Survey (MzLS; Prop. ID #2016A-0453; PI: Arjun Dey). DECaLS, BASS and MzLS together include data obtained, respectively, at the Blanco telescope, Cerro Tololo Inter-American Observatory, NSF’s NOIRLab; the Bok telescope, Steward Observatory, University of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab. The Legacy Surveys project is honored to be permitted to conduct astronomical research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. NOIRLab is operated by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo, Financiadora de Estudos e Projetos, Fundacao Carlos Chagas Filho de Amparo a Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Cientifico e Tecnologico and the Ministerio da Ciencia, Tecnologia e Inovacao, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energeticas, Medioambientales y Tecnologicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenossische Technische Hochschule (ETH) Zurich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciencies de l’Espai (IEEC/CSIC), the Institut de Fisica d’Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig Maximilians Universitat Munchen and the associated Excellence Cluster Universe, the University of Michigan, NSF’s NOIRLab, the University of Nottingham, the Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University. The Legacy Survey team makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer (NEOWISE), which is a project of the Jet Propulsion Laboratory/California Institute of Technology. NEOWISE is funded by the National Aeronautics and Space Administration.

References

Appendix A: Selection of INS candidates

To search for INS candidates, we started with a preliminary screening to remove all spurious sources arising from soft strips that appear in some parts of CCD4 in MOS1. We then selected sources with positions in the HR₂-HR₃ hardness ratio diagram consistent with those of the M7 group of X-ray thermally emitting INSs (see, e.g. Fig. 1 of Pires et al. 2009 and similar hardness ratio diagrams in Rigoselli et al. 2022b; we refer to the caption of Table 1, for a definition of hardness ratios and standard XMM-Newton energy bands). Practically, this implies selecting sources with a very soft spectral slope in the 0.5 − 2 keV range and consistent with no detected flux above 2 keV, namely: HR₂ − σ_HR2 < −0.3 and HR₃ − σ_HR3 < −0.99. Since HR₁ mostly depends on interstellar absorption, we did not put any constraint on its value. In order to exclude optically bright classes of soft X-ray emitters such as active coronae or cataclysmic variables, we removed XMM-Newton sources having a positional match with a probability higher than 50% with a selection of archival catalogues (Pineau et al. 2017). For this purpose we used the results of the statistical cross-matching process between 4XMM-DR9 and GALEX GR6+7, XMM-OM-SUSS4.1, Gaia DR2, APASS, SDSSDR12, Pan-STARRS DR1, SkyMapper, 2MASS and AllWISE (see Webb et al. 2020, for a description). Visual screening allowed us to discard high proper motion active coronae which cannot match the X-ray source position due to epoch difference. Finally, we looked for possible identifications in the SIMBAD⁵ astronomical database.

Appendix B: Inter-epoch variability in eRASS

The best-fit parameters resulting from the analysis of the eROSITA epochs, when fitted independently of EPIC, have shown a strong dependence on the choice of TM detectors, patterns, and binning. For example, adopting a simple blackbody model, the best-fit temperature of the target derived from TM8 is 30% softer than when TM 5 and 7 are included in the analysis. More consistent results between the eROSITA detectors were obtained when N_H is fixed to the best-fit EPIC value. In particular, the fit of a constant blackbody model on the eROSITA data sets has a high null-hypothesis probability, 80% for 19 degrees of freedom, and parameters in overall agreement with EPIC: $kT={67}_{-6}^{+7}$ eV, $R_{\mathrm{em}}=7_{-2}^{+3}$ km and $f_{\mathrm{X}}=(1.{15}_{-0.19}^{+0.15})\times {10}^{-13}$ erg s⁻¹ cm⁻².

For model (4) in Table 4 we only kept the column density towards the target and the parameters of the blackbody within a given epoch (that is, between the EPIC cameras) coupled during fitting. When the requirement of constancy is lifted from the individual epochs, the parameters of the target as derived from eRASS in model (4) show a much wider spread around the median in comparison to the EPIC data, with variations in temperature and flux of up to 13% and 60%, respectively (Table 5 and Fig. 4). The chi-square values assuming constancy are 3.6, 13.2, and 8.2 (5 d.o.f.), respectively for kT, R_em and f_X.

We investigated whether similar variations are observed in the survey data of other thermally emitting INSs. We retrieved the c020 pipeline processed sky tiles centred on the positions of RX J0420.0−5022, RX J0720.4−3125, RX J0806.4−4123, RX J1308.6+2127, RX J1856.5−3754, and 2XMM J104608.7−594306 (Table B.1), and analysed them consistently with the procedure described here for J0221. The 0.2 − 1 keV count rates of the sources from PSF fitting as a function of eRASS epoch can be seen in the diagram of Fig. B.1, with 1σ errors. A chi-square test for constancy at p-value 0.1 is formally rejected for all sources except RX J1856.5−3754 and 2XMM J104608.7−594306 (χ² values within 1.9 and 6.8 above the critical value of ${\chi }_{\text{cr}}^2\sim 0.58$ for 3 d.o.f.). All rates are consistent within 2σ.

Fig. B.1. Count rates of thermally emitting INSs in eRASS1–4 (TM8, 0.2 − 1 keV). A chi-square test for constancy is formally rejected for 4XMM J022141.5−735632, RX J0420.0−5022, RX J0720.4−3125, RX J0806.4−4123, and RX J1308.6+2127. All rates are consistent within 2σ.

To fit the individual eRASS1–4 spectra of the six INSs we assumed simple absorbed blackbody models. The fit results are in general agreement with those previously reported in the literature. We found up to 2σ variations in flux and deviations of up to 23% in kT, suggesting that the level of intra-eRASS variability observed for the INS candidate J0221 is typical and probably a result of low-count statistics.

Appendix C: Variable elemental abundances

In Table C.1, we show the results of two additional spectral fits with the variable absorption model tbvarabs in XSPEC, which were not included in Table 4 for conciseness. Model (A) is a single-temperature constant blackbody with enhanced oxygen abundance with respect to solar, $Z_{\mathrm{O}}=(5.4_{-1.9}^{+2.4})\times Z_{\odot }$; it consists of a significant improvement with respect to model (1) in Table 4 (F-test probability 0.003). Elemental abundances other than oxygen were either insensitive or unconstrained by the fit. Similarly to the improvement found for the multi-component models (2) and (3), the best-fit parameters of the blackbody component do not significantly (> 2σ) differ from those of model (1).

In model (B), we tested the possibility that the source is extragalactic, motivated by its projected location in the Magellanic Bridge and N_H in excess of the Galactic value from Dickey & Lockman (1990), $N_{\mathrm{H}}^{\mathrm{DL}}=3.8\times {10}^{20}$ cm⁻². As discussed in Sect. 4, the hydrogen column density integrated in the line-of-sight includes the contribution of both the Milky Way and the Magellanic Clouds; ideally, these should be uncoupled to get a reliable estimate of the foreground N_H in our Galaxy (e.g. Sasaki et al. 2022). Assuming that the source is located at the distance of the Small Magellanic Cloud (SMC), d_SMC = 60 kpc, we thus tested a model with two absorption components. We set $N_{\mathrm{H}}^{\mathrm{Gal}}\equiv N_{\mathrm{H}}^{\mathrm{DL}}$ whereas $N_{\mathrm{H}}^{\mathrm{SMC}}$ was free to vary. For the latter, the elemental abundances of Anders & Grevesse (1989) were initially fixed to a sub-solar value of 0.2, more typical of that of the SMC (Russell & Dopita 1992), then allowed to vary to search for an improvement in the fit statistic. The best solution, with $N_{\mathrm{H}}^{\mathrm{SMC}}=3.4_{-2.0}^{+2.2}\times {10}^{20}$ cm⁻² and kT = 67.2 ± 2.9 eV, is obtained when oxygen is overabundant with respect to solar, $Z_{\mathrm{O}}=(4_{-2}^{+8})\times Z_{\odot }$; we note that this is essentially the same fit result as that of model (A). In Table C.1 we list the corresponding emission radius of the model normalised to d_SMC, $R_{\mathrm{em}}^{\mathrm{SMC}}={340}_{-50}^{+70}$ km.

All Tables

All Figures

	Fig. 1. Smoothed false-colour (top) and source detection (bottom) images of the two XMM-Newton observations covering the position of J0221. The red band covers the energies from 0.2−1.0 keV, green from 1.0−2.0 keV, and blue from 2.0−12.0 keV. The position of the target and of other sources detected in the field-of-view are marked by white and light blue circles, respectively.
In the text

	Fig. 2. eROSITA extraction regions (eRASS:4, TM0, 0.2 − 1 keV). Top panel: background annulus (magenta) with inner and outer radii of sizes 2.5′ and 15′, respectively, excluding “contaminants” (cyan circles with a red bar across). Bottom panel: source circular region (yellow) of size 82″.
In the text

	Fig. 3. Best-fit double-blackbody model and residuals (fit 2 in Table 4). The two EPIC pn and four MOS data sets are colour-coded in dark grey and blue, respectively. The eROSITA spectrum (merged here for display purposes) is shown in pink.
In the text

Fig. 4. Long-term spectral and flux variation of the X-ray source. Left panel: blackbody temperature and emission radius as a function of MJD (“variable” model 4; see Tables 4–5). Right panel: long-term evolution of the 0.2 − 2 keV X-ray flux of the target, including upper limits and previous detection by other X-ray missions. The time interval extends back to the ROSAT All-Sky Survey (RASS) and pointed (PSPC) era and include data points from Swift XRT and XMM-Newton slew observations (see the text, for details). In all plots the purple horizontal shaded areas show the 1σ median absolute deviation of the parameters.

In the text

	Fig. 5. Photometric-weighted multi-filter stack of Legacy Survey DR10 g, r, z and DECam u, i images centred on the sky position of J0221. The blank (mG > 26.45) 3σ confidence level error circle of the X-ray source is displayed in the inset in an inverted black and white colour map.
In the text

	Fig. 6. Thermal luminosity vs. emission radius for cooling neutron stars. The arrow shows the distance-luminosity range of the target assuming the same emission radii observed for the M7 INSs (see the text, for details).
In the text

	Fig. B.1. Count rates of thermally emitting INSs in eRASS1–4 (TM8, 0.2 − 1 keV). A chi-square test for constancy is formally rejected for 4XMM J022141.5−735632, RX J0420.0−5022, RX J0720.4−3125, RX J0806.4−4123, and RX J1308.6+2127. All rates are consistent within 2σ.
In the text