© The Authors 2025

1. Introduction

In recent years, there has been renewed interest in identifying X-ray emitting isolated neutron stars (XINSs) that are not detected by conventional radio and gamma-ray pulsar surveys (Rigoselli et al. 2022; Pires et al. 2022; Demasi et al. 2024; Kurpas et al. 2023, 2024b). While the detectability of old and distant radio pulsars may be affected by several selection biases, including beaming, nulling, and propagation effects in the interstellar medium (ISM), gamma-ray surveys such as those conducted by the Fermi Gamma-ray Space Telescope (Atwood et al. 2009) can miss pulsars if their gamma-ray emission is weak or if the pulsars are located in regions of strong gamma-ray background. X-ray observations, though affected by absorption at soft energies by the ISM, are sensitive to both thermal emission from the hot neutron star surface and the magnetospheric radiation powered by neutron star rotation. In particular, searches for thermally emitting XINSs provide a detection channel that is independent of the pulsar’s spin and viewing geometry.

Wide-field X-ray survey missions, such as ROSAT in the 1990s (Truemper 1982) and eROSITA on board the recently launched Spectrum Roentgen Gamma (SRG) Observatory (Predehl et al. 2021), greatly enhance our understanding of neutron star demographics, their physical properties, and their role in the broader astrophysical context. The X-ray source eRASSU J065715.3+260428 (hereafter dubbed J0657; Kurpas et al. 2023, 2024b) is one of the first candidates identified in a search for XINS in the western Galactic hemisphere of the eROSITA All-Sky Survey (eRASS; Merloni et al. 2024). Its X-ray spectrum, best described by a single blackbody component with kT ∼ 110 eV, is predominantly thermal and seemingly constant. Deep optical follow-up imaging with the Large Binocular Telescope (Hill et al. 2012) in the r band has excluded the presence of counterparts brighter than 25.5 mag, suggesting a compact nature (X-ray-to-optical flux ratio f_X/f_opt ≳ 860).

Additional observations are needed to fully characterise the evolutionary state of the source and place it among the families of Galactic isolated neutron stars (INSs; e.g. Pires et al. 2017; Khokhryakova et al. 2021; Popov 2023, for an overview). For example, the measurement of the INS rotation and spin evolution provides estimates of the magnetic field strength and characteristic age under the usual assumptions of a dipolar magnetic field configuration and loss of rotational energy solely due to magnetic braking in vacuum (Ostriker & Gunn 1969). Spectral characterisation is equally important to constrain the presence of possible absorption features and non-thermal components. However, the shallow survey exposures of eROSITA (typically, ∼150 s per sky visit for this source), and its comparably low effective area at energies above 1 keV, prevent detailed timing and spectral analysis.

We report in this work the results of a dedicated multi-wavelength follow-up campaign targeting the XINS J0657. At X-ray energies, J0657 was observed with the Neutron star Interior Composition Explorer (NICER; Gendreau et al. 2016) and XMM-Newton (Jansen et al. 2001). In the optical, we obtained deep imaging with the Focal Reducer/low dispersion Spectrograph 2 (FORS2) instrument at the European Southern Observatory Very Large Telescope (ESO-VLT; Appenzeller et al. 1998). In the radio regime, we performed a DDT observation with the Five-hundred-meter Aperture Spherical radio Telescope (FAST; Nan et al. 2011), which was carried out simultaneously with the XMM-Newton pointing. We complemented the analysis with gamma-ray observations from the Fermi Large Area Telescope. The XINS candidate J0657 is the second eROSITA selected X-ray source to be the target of such a multi-wavelength effort. We previously reported the results of follow-up observations of eRASSU J131716.9–402647, a long spin-period X-ray source displaying properties that closely resemble those of highly magnetised XINSs (Kurpas et al. 2024a). These first promising results confirm eROSITA’s ability to unveil elusive XINSs beyond the solar neighbourhood.

The outline of the paper is as follows: we begin with a description of the observations and data reduction in Sect. 2. The results of the X-ray timing and spectral analysis, along with the search for optical, radio, and gamma-ray counterparts, are presented in Sect. 3. In Sect. 4, we discuss the implications for the nature of J0657, particularly in the context of the known population of Galactic INSs. Finally, our conclusions are provided in Sect. 5.

2. Observations and data reduction

2.1. NICER

Between November 9 and 25, 2023, NICER observed the XINS candidate J0657 continuously for approximately 210 ks. We used the NICERDAS tools version 2024-02-09_V012A, distributed as part of the HEASoft software release version 6.33.2, to process the observations. We applied the nicerl2 task to extract the event lists and remove flagged events. We only considered intervals of orbital night-time to mitigate the impact of the optical leak affecting the scientific performance of NICER¹. As J0657 is a fairly faint and soft X-ray source, we adopted conservative screening criteria according to the NICER analysis guidelines². The cleaned observations were merged by the niobsmerge task. Periods of remaining high flaring background were removed through a 3σ-clipping algorithm from the mean count rate in the 0.3 − 12 keV energy band. The total net exposure after applying all screening criteria is 116 ks and the total events ∼8 × 10⁴ (0.3 − 2 keV; Table 1). Spectra and light curves were extracted with the nicerl3-spect and nicerl3-lc pipelines. We estimated the background with the SCORPEON³ model. For the timing analysis, the barycorr task, JPLEPH.405 ephemeris, and the improved X-ray sky position of the target derived from the XMM-Newton EPIC observation (Sect. 2.2) were used to correct the times-of-arrival (ToAs) to the barycentre of the Solar System.

2.2. XMM-Newton

As part of large programme 092128, the XMM-Newton Observatory targeted the XINS J0657 on March 14/15, 2024, for ∼67 ks (Table 1). The EPIC pn detector (Strüder et al. 2001) was operated in full-frame (FF) mode, whereas the small-window (SW) mode was selected for EPIC MOS (Turner et al. 2001). Both cameras were equipped with the thin filter for its enhanced response at soft X-ray energies. We used the XMM-Newton Science Analysis Software (version: 21.0.0) to process the observations. We found the observation to be significantly affected by periods of high background (about 30% of the elapsed time). Similarly to NICER (Sect. 2.1), we applied a 3σ-clipping algorithm to remove time intervals of flaring background. The total net exposure is ∼45 − 50 ks for the individual EPIC cameras, amounting to a total 1.0 × 10⁴ counts (0.2 − 2 keV; Table 1).

We conducted source detection with the edetect_stack task across the five standard XMM-Newton bands and three EPIC cameras (Traulsen et al. 2020). The astrometry was refined with the SAS task eposcorr, considering 27 potential matches of common sources from the optical Guide Star Catalogue (version: 2.4.2; Lasker et al. 2008) that are located within 15′ of the nominal pointing coordinates of the EPIC observation. The resulting X-ray sky position is given in equatorial and Galactic coordinates in Table 2. Based on the signal-to-noise (S/N) ratio, optimal source and background regions were defined with the SAS task eregionanalyse in the 0.2 − 12 keV energy band.

We followed the general SAS guidelines to extract light curves and spectra. Single and double pattern events were considered for pn (PATTERN < = 4) while all valid patterns were included for the MOS detectors (PATTERN < = 12). The extracted spectra were grouped using a maximum oversampling factor of three and at least 25 counts per energy bin. The barycentric correction was performed with the SAS barycen task, DE405 ephemeris, and the updated X-ray source position of the target listed in Table 2.

2.3. ESO-VLT

The field of J0657 was observed with the FORS2 instrument at the ESO-VLT in December 2023 for 1 hour. In total, nine images in the R_SPECIAL band were taken, each with 400 s exposure time. We used the EsoReflex environment with the FORS workflow for imaging data (version 5.6.2) to conduct the basic data reduction (Freudling et al. 2013). The resulting images were subsequently astrometrically calibrated using the astrometry.net code (Lang et al. 2010). As several bright stars were saturated in the FORS images, an individual index file based on Legacy Survey DR10 (Dey et al. 2019) positions was adopted for this task. To facilitate stacking, the individual exposures were reprojected with the astropy reproject package (Astropy Collaboration 2013, 2018, 2022). Finally, we removed cosmic rays with the L.A.Cosmic algorithm (van Dokkum 2001; McCully et al. 2018) implemented in the ccdproc package (Craig et al. 2017).

2.4. FAST

To detect potential periodic signals in the radio band for J0657, FAST was used with the L-band 19-beam receiver (Jiang et al. 2020) for the observations. The system temperature is approximately 24 K and the central beam size at 1.4 GHz is 2.9′. Given a pointing error of less than 8″ and the target J0657 having a position error of 0.2″ (see Table 2), it should be well within the beam and very close to the centre. Observations were conducted on March 14, 2024, for a duration of 4.86 hours. Data were sampled with 8-bit precision for four polarisations, with a sampling time of 98.304 μs and divided into 1024 channels covering 1–1.5 GHz (with a channel width of 488 kHz). The data were recorded in PSRFITS format (Hotan et al. 2004) for subsequent analysis.

3. Results

3.1. X-ray timing analysis

We searched for periodic modulations in the NICER and XMM-Newton EPIC pn datasets applying a Lomb-Scargle algorithm (Lomb 1976; Scargle 1982). We adopted time bins of 1 ms for NICER and 90 ms for the XMM-Newton light curves. We found both time series to contain a significant (> 4σ − 8σ) periodic signal at ∼261 ms (see Fig. 1 and Table 3). While the first harmonic is detected at twice the fundamental frequency in the NICER dataset, it cannot be detected in the XMM-Newton data. Only a significant signal at 137 ms is identified, that arises due to a beating between the 90 ms binning and the 261 ms modulation of the pulsar. We searched for shorter period modulations in the NICER dataset (down to 1 ms), but found no additional significant signal. The EPIC MOS observations were not included in the search, due to their insufficient time resolution⁴.

Fig. 1. Lomb-Scargle periodogram for the NICER observation. The horizontal green dashed lines show the power corresponding to significance levels of 1σ to 8σ.

The pulse profile in Fig. 2, folded at the fundamental peak at 261 ms, shows a single pulse and a saddle between phases 0.2 and 0.6 of the pulsar rotation. Folding the data at the first harmonics instead, leaves only a single broad peak for roughly half the pulsar rotation that can generally be well modelled with a single sinusoidal function (reduced chi-squared χ_ν² ∼ 1 for ν = 12 degrees of freedom in the NICER dataset). Since multiple harmonics are necessary to accurately model the pulse profile in Fig. 2, we interpret the 261 ms modulation as representing the true rotational period of J0657. This is because it includes more distinct and well defined features.

Fig. 2. X-ray pulse profiles from NICER (top) and XMM-Newton (top), folded at the 261 ms modulation (see Table 3 and the text, for details). The two phase bins that were defined for the phase-resolved analysis of the XMM-Newton data are marked by the light- and dark-grey regions in the lower panel.

With the goal to refine the period estimation for a phase-coherent analysis, we applied the Gregory & Loredo (1996) method (see e.g. Kurpas et al. 2024a, for a description). Specifically, we adopted m_max = 12 and frequency steps (Δν_GL) of 1 μHz and 20 μHz for NICER and XMM-Newton, respectively. The results are shown in Table 3. The spin period and pulsed fraction between XMM-Newton and NICER agree within 2σ; we note that the deviations between the datasets may be attributed to the stronger NICER background and the higher timing accuracy in comparison to EPIC pn.

The mid-observation times of NICER and XMM-Newton are roughly separated by 120 days. The linear timing solution from NICER (arbitrarily centred at the phase of the maximum of the pulse profile, see Table 3), indicates that both observations can be connected in phase with a cycle counting error of 0.25. However, folding the XMM-Newton events with the NICER timing solution reveals a phase shift of 0.4 between the maxima of the pulse profiles, indicating the presence of a measurable spin-down between the two datasets.

To find a best-fitting phase coherent timing solution, we defined the residuals between observed and tested cycle counting as

$$\begin{array}{c}\hfill R=\sum _{i=1}^N{\left(\frac{C_{\mathrm{obs},\mathrm{i}}-C(t_i,\nu ,\dot{\nu },t_0,C_0)}{\sigma (C_{\mathrm{obs},\mathrm{i}})}\right)}^2,\end{array}$$(1)

where the summation is over the N defined ToAs of the events. In Eq. (1), we adopted a quadratic timing solution of the form:

$$\begin{array}{c}\hfill C(t,\nu ,\dot{\nu },t_0,C_0)=C_0+\nu (t-t_0)+\frac{1}{2}\dot{\nu }{(t-t_0)}^2,\end{array}$$

where ν and $\dot{\nu }$ is the tested spin and spin-down, t₀ is the reference time, and C₀ is the phase shift.

We split the NICER exposure into sections of at least 13 000 photons, resulting in a series of six ToAs where the XINS pulsation could be reliably detected. To determine the phase of each pulse, we fitted the pulse profile (folded with the tested timing solution) with a sinusoidal function including the fundamental and first harmonic. The fit parameters were then used to estimate the phase of the maximum of the pulse profile and, consequently, the ToA that is closest to midpoint of each time section. We then searched for the timing solution that minimised the cycle counting residuals for our data.

We explored the likelihood landscape of possible timing solutions through the nested sampling algorithm implemented in the UltraNest package (Buchner 2021). We assumed linear priors for the period and phase shift, whereas a logarithmic prior was used for the spin-down. The reference time was kept near the mid-observation of the NICER dataset and near the time of phase zero at t₀ = 60265.4441739. In the analysis, we evolved 800 sample points until the remaining unexplored fraction of the evidence integral amounted to less than 1%.

The resulting posterior probability distribution is shown in Fig. 3. We obtained from the distribution the best solution P = 261.085400(4) ms, $\dot{P}=6_{-4}^{+11}\times {10}^{-15}$ s s⁻¹, and $C_0=0.{001}_{-0.012}^{+0.013}$ (we give the median value and confidence range, containing 68% of the sample points). We note, however, that this timing solution still allows for ambiguity in the cycle counting. While the solution constrains the possible range of values for P and Ṗ, there are still multiple points within the confidence region space that can minimise the cycle counting residuals, but may deviate by multiple cycles from each other. Additional observations covering a wider sampling in time will be needed to correctly identify the best-fit solution from these low residual points.

Fig. 3. Corner plot depicting the parameter distribution of the coherent timing search. For presentation purposes, we subtracted the period P0 = 261.085 ms from the period values of the distribution. Vertical lines mark the 15.8%, median and 84.1% percentiles.

3.2. X-ray spectral analysis

3.2.1. Phase-integrated analysis

The spectral analysis of the XMM-Newton and NICER datasets was performed with XSPEC (version: 12.14.0h, Arnaud 1996). All models were coupled with a tbabs component to account for interstellar absorption (adopting elemental abundances from Wilms et al. 2000). While analysing the spectra from multiple instruments simultaneously (XTI, EPIC pn, and MOS1/2), we included an additional constant factor in the fit to address cross-calibration uncertainties. Chi-squared statistics were employed to fit the spectral data.

The results are presented in Table 4. We initiated the spectral analysis using the EPIC pn dataset for simplicity. Initially, a single-component absorbed blackbody (BB) model did not provide a satisfactory fit (χ_ν²(ν) = 5.29(42); see fit I in Table 4). Including a second BB component improved the fit residuals compared to fit I (cf. fit II in Table 4). Alternatively, a simple absorbed power-law model (PL; fit IV in Table 4) closely matched the data. However, in this scenario, the best-fit column density N_H exceeds the Galactic value along the line of sight, N_H, gal = 8.27 × 10²⁰ cm⁻² (HI4PI Collaboration 2016), and the photon index α of the power-law model is much steeper than typically observed in rotation-powered pulsars (RPPs; Becker et al. 2009).

Absorption features have been observed in the spectra of several XINSs (e.g. Tiengo et al. 2013; Schwope et al. 2022; Pires et al. 2023). Similarly, for J0657, we found that including a Gaussian absorption component (GABS) at low energies (∼0.3 keV) significantly enhances the fit quality (models III and VI in Table 4), although there exists a strong degeneracy between the line width σ and N_H for model III. No improvement was observed when σ was held fixed during fitting. For model VI, the power-law spectral slope remains significantly steeper than what is typically observed in the known INS population. We present the folded spectrum in Fig. 4 along with the best-fit 2BBGABS model. For comparison, the fit residuals using single or double BB components are also shown. For an assumed distance of 1 kpc, we derive the size of the emission region for the dominant colder BB component in models I, II, and III. With values ranging from 1 to 4 km the emission region size is below the canonical neutron star radius of ∼12 km, but still consistent with phenomenological models used to describe the spectra of the purely thermally emitting XINSs (e.g. the ‘magnificent seven’ described by e.g. Turolla 2009; Potekhin et al. 2020).

Fig. 4. EPIC pn spectrum of J0657 along with the best-fit single blackbody (red, I), two blackbody (blue, II) or two blackbody with a line (green, III) models. See Table 4 for the fit parameters.

The fit of a composite model combining thermal and non-thermal components reveals a steep PL and a hot BB component (fit result VII in Table 4), contrasting with middle-aged RPPs where soft thermal components dominate (e.g. the ‘three musketeers’ described by De Luca et al. 2005). In these sources, X-ray emission is often interpreted as a blend of thermal radiation from the cooling surface and hot polar caps. Additionally, fainter power-law tails extending to higher energies are attributed to magnetospheric activity of the pulsar.

We calculated upper limits for the detection of comparable non-thermal components using a 2BBGABSPL model. We held constant the parameters of the BB and GABS components at the optimal values from model III, while assuming a PL slope of α = 2. The normalisation of the power-law component was allowed to vary freely. This analysis yielded a 3σ upper limit of 2.1 × 10⁻¹⁴ erg s⁻¹ cm⁻², equivalent to about (2.1 ± 1.8)% of the source’s unabsorbed flux of model III.

Qualitatively similar outcomes were observed using neutron star atmosphere (NSA; Pavlov et al. 1995; Zavlin et al. 1996) models (fit results VIII, IX, and X in Table 4). A single NSA component alone fails to adequately fit the data, regardless of the magnetic field intensity. However, composite models incorporating the 0.3 keV absorption feature, alongside either a hot BB component with kT ∼ 250 eV or a steep PL with α ∼ 4.6, achieve statistically satisfactory fits. The temperature and normalisation of the atmosphere are notably affected by the inclusion of additional components: in the NSAGABSPL scenario, we observe a hot atmosphere and a substantial distance, while in the NSAGABSBB scenario, a colder atmosphere and a closer source are inferred.

To assess the significance of the absorption feature, we calculated both the false-negative (where an existing feature is missed or incorrectly identified) and false-positive (where a non-existing feature is identified) rates. For the false-negative rate, we utilised the XSPECfakeit command to simulate 1000 spectra based on the best-fit result from Table 4 (III), employing a 2BBGABS model. We then performed fits using both a 2BB model and the 2BBGABS model to determine how often the simulated absorption line would be undetected. We found a low false-negative rate of only 6.2%, indicating that the spectral parameters were generally well recovered and confirming that the EPIC pn observation accurately captures the source’s spectral properties.

For the false-positive rate, we simulated 1000 spectra using the best-fit solution from Table 4 (II), employing a 2BB model. These simulated spectra were then fitted using both a 2BB model and a 2BBGABS model. We found that in only 2.7% of all simulations did the addition of the absorption line improve the fit statistic by more than the observed Δχ_ν² > 0.14. This low false-positive rate suggests a high significance that the observed feature is real.

The data from additional X-ray instruments, EPIC MOS1/2 and NICER XTI, were next included in the analysis. Simultaneous fits with all EPIC instruments aligned with previous results (fit XI in Table 4). However, incorporating NICER data was challenging, as the modelling requires correct handling of background contamination, especially if bright field sources are present in the vicinity of the target. In particular, we identified in the EPIC images an X-ray source of similar brightness to the XINS, 2.10(5)×10⁻¹³ erg s⁻¹ cm⁻² (0.2 − 10 keV), 2′ north-west of its position. Its EPIC pn spectrum is best-fit by an absorbed PL with $N_{\text{H}}=7.{20}_{-1.0}^{+1.1}\times {10}^{20}$ cm⁻² and α = 1.98(6). To address this contamination and the large background at soft energies, a power-law component was added to the spectral model, and only NICER photons above 500 eV were included. The final simultaneous spectral fit aligned well with EPIC pn results (XII), except for some flux discrepancies likely due to background contributions or calibration inaccuracies.

3.2.2. Phase-resolved analysis

Achieving a detailed phase-resolved analysis with the current dataset is challenging. For EPIC pn, the frame time in FF mode is a significant fraction of the spin period of J0657, leading to some event confusion even with only three phase bins. Although NICER data could, in principle, provide a more detailed phase-resolved analysis, the increased background and contamination led us to focus on the EPIC pn data.

Within these limitations, we defined two phase intervals aimed at distinguishing the pulsed (ϕ_pulse = 0.85 − 0.2) and non-pulsed (ϕ_{off − pulse} = 0.2 − 0.85) emission of J0657 (see the grey background in the lower panel in Fig. 2). For each phase bin, we extracted the source spectrum. We did not correct the exposure for the camera dead time⁵.

We conducted simultaneous fits of the pulsed and non-pulsed spectra, fixing only the N_H value and σ width of the line and leaving other model components free. The results were similar to those for the phase-integrated spectrum, with PL components converging to unreasonably large slopes and requiring at least two BB components to fit the spectrum. For the double blackbody models in Table 5, the additional GABS absorption component was not necessary, although it improved the fit slightly and led to smaller N_H values. The model parameters showed a maximum 2σ variation between phase bins, which is not significant. The absolute temperatures of both BB components seemed slightly higher during the pulse maximum. Additional observations are needed to confirm any phase-variable trends.

3.3. Searches for pulsed radio emission

The pulsar search, which included both single pulse search and periodic signal detection, was conducted using PRESTO V4.0⁶ (Ransom 2011). The Radio Frequency Interference (RFI) was firstly masked by the rfifind routine, with the search for RFI in the frequency domain conducted over a time span of 12.9 seconds. Then, the data were de-dispersed by the prepsubband routine, which were corrected for the dispersion effects based on the pulsar’s position and distance (measured in kpc). We used the predicted dispersion measure (DM) from the model by Yao et al. (2017, YMW16). According to this model, the predicted DM is < 200 pc cm⁻³. To ensure we captured any potential signals, we set the upper limit of the DM range to twice the predicted value, which is ∼400 pc cm⁻³, and the DM was incremented in steps of 0.1 pc cm⁻³ across the entire range from zero.

The de-dispersed time series were used for further single pulse searches, or transformed into the frequency domain via Fast Fourier Transform (FFT) for periodic signal detection. For the single pulse search, the PRESTO routine single_pulse_search.py was applied. The accelsearch routine was used to search for periodic signals. Although J0657 does not appear to be a binary pulsar, we set the maximum allowable frequency shift in the frequency domain (measured in Fourier bins) to 50 using the ‘-zmax’ option. The candidates were produced for further identification and confirmation by both the ACCEL_sift.py routine from PRESTO and manual selection after determining the possible spinning periods. None of the candidates were identified as either J0657 or any other pulsar.

The sensitivity of our search can be estimated using the radiometer equation (e.g. Lynch et al. 2011):

$$\begin{array}{c}\hfill S_{\min }=\beta \frac{(S/N_{\min })T_{\mathrm{sys}}}{G\sqrt{n_{\mathrm{p}}{t}_{\mathrm{int}}\mathrm{\Delta }f}}\sqrt{\frac{W}{P-W}},\end{array}$$

where β represents the sampling efficiency and is equal to 1 for our 8-bit recording system. In this equation, W denotes the width of the pulsar profile; we used 10% of the pulse width for our calculations. The minimum signal-to-noise ratio, (S/N)_min, for our search is set to 10. The system temperature, T_sys, is 24 K, and the antenna gain, G, is 16 K Jy⁻¹. The number of polarisations, n_p, is 2. The integration time t_int is measured in seconds, and the bandwidth Δf is 300 MHz. This bandwidth accounts for the fact that approximately 25% of the data were masked as RFI and thus excluded from further processing. Consequently, the radio L-band flux of J0657 should be lower than 1.4 μJy (10σ).

3.4. Optical limit

The ESO-VLT observation of J0657 in the R_SPECIAL band provides a deeper optical limit than previously achieved with LBT (Kurpas et al. 2023). In Fig. 5, we mark the positions of nearby field sources detected using SEXTRACTOR (Bertin & Arnouts 1996), as well as the position of the target determined from eROSITA (blue circle) and XMM-Newton (green circle; Table 2). Clearly, there is no sign of a source in this blank region: the closest optical object is at a > 8σ separation from the well-determined EPIC pn position of the X-ray source.

Fig. 5. FORS2 R-band image of the field of J0657. The eROSITA (Kurpas et al. 2023, 2024a) and XMM-Newton X-ray sky positions are marked by the blue and green circles, respectively (radius marks the 1σ positional accuracy). Brown circles (with arbitrary radius of 1″) mark field sources identified from an SExtractor run.

The 5σ detection limit of the FORS2 image is 27.32 mag (Table 6). Using the relation presented in Sect. 3.4 of Kurpas et al. (2024b) and the best-fit 2BBGABS model flux of Model III in Table 4, we determined the X-ray-to-optical flux ratio limit to be above 5260, an improvement of about one order of magnitude compared to the previous LBT imaging.

3.5. Gamma-ray analysis

No counterpart is present within ∼1.1° of the position of the X-ray source in the fourth Fermi Source Catalogue Data Release 4 (4FGL-DR4; Abdollahi et al. 2022; Ballet et al. 2023). The closest entry is the blazar candidate 4FGL J0658.2+2709, with an energy flux of 1.8(4)×10⁻¹² erg s⁻¹ cm⁻² (1σ) in the 0.1 − 100 GeV energy range.

To search for gamma-ray emission from J0657, we retrieved Pass 8 Fermi-LAT photon data accumulated between 2008 August 4 and 2024 June 11 within a region of interest (ROI) of 15° from the XINS position. The analysis was carried out with fermitools⁷ following the guidelines from the Fermi Science Support Center. The photon files were filtered using gtselect within 0.1 − 500 GeV and zenith angles smaller than 90°. We only considered ‘source’ events with evclass=128 and evtype=3. Good-time intervals were applied with gtmktime and the filter expression (DATA_QUAL==1) && (LAT_CONFIG==1), ensuring that only data with the highest quality were used. Additionally, events registered during time intervals identified in the GTI extension of the 4FGL-DR4 catalogue as being associated with solar flares and gamma-ray bursts were removed from the analysis.

We carried out a binned likelihood analysis with gtlike to determine whether the XINS is detected as a point source (Abdo et al. 2009). For this purpose, we adopted the recommended Galactic diffuse emission model, gll_iem_v07.fits, and isotropic spectral template, iso_P8R3_SOURCE_V3_v1.txt, and created a source model XML file with the locations, spatial and spectral shapes of the 4FGL sources within 25° from the XINS, which may contribute photons to the analysis at different energies. Count and source maps as well as the mission exposure and livetime were generated using gtbin, gtsrcmaps, gtexpcube2 and gtltcube.

The normalisation parameters of all bright and variable sources within 10° of the X-ray position, as well as the spectral parameters of the diffuse components, were left free during fitting; the parameters of the remaining 4FGL sources were kept fixed at their catalogued values. Altogether, the contribution of 277 Fermi sources, including four extended, was taken into account in the computation. Finally, we adopted a two-step minimisation procedure with gtlike, which is recommended for faint sources. Initially, we fitted the data with a tolerant convergence criterion using the MINUIT optimiser. This was followed by a more precise parameter estimation with NEWMINUIT over the initial results, using a strict fit convergence threshold.

Before introducing a putative source at the position of the XINS, we first performed fits to replicate the flux of sources catalogued in the 4FGL-DR4, which is based on a slightly shorter time baseline, and to achieve flat residuals across all energies. As an initial approximation, we then modelled the spectrum of the tentative gamma-ray counterpart of the INS using a simple power-law over the full sixteen-year time baseline of the analysis.

Consistent with the 4FGL-DR4 results, the likelihood analysis yields a low test statistic (TS) value of 6 at the position of the XINS, confirming the absence of a gamma-ray counterpart. Using the best-fit likelihood model, we derived upper limits for the pulsar’s photon flux with the UpperLimits Python module in fermitools, allowing both the power-law normalisation and photon index to vary freely. The resulting 95% confidence upper limit, 7 × 10⁻⁹ photons s⁻¹ cm⁻², corresponds to an energy flux of 6.8 × 10⁻¹³ erg s⁻¹ cm⁻² in the 0.1 − 100 GeV energy range.

Despite the lack of a spatial detection, we searched for pulsed gamma-ray emission using the X-ray timing solution from the XMM-Newton and NICER campaigns (see Sect. 3.1). We used an aperture radius of 5° for events above 100 MeV. The photons were barycentred using the source coordinates listed in Table 2 and the JPL DE405 solar system ephemeris. The short time baseline of the X-ray ephemeris (∼126 days) resulted in too few gamma-ray photons (∼5000) for a meaningful timing analysis. Therefore, we extended the search to cover the entire time baseline of the analysis. We adopted conservative ±3σ ranges in ν and $\dot{\nu }$ to account for all plausible values. To assess the significance of a potential signal, we used an implementation of the H-test optimised for detecting faint gamma-ray pulsars in Fermi-LAT data⁸. This approach employs the ‘simple weights’ approximation described by Bruel (2019), Smith et al. (2019). The searches yielded no significant modulation.

4. Discussion

In this work, we report the results of a dedicated multi-wavelength follow-up campaign targeting the XINS candidate J0657. The deep ESO-VLTR-band observation confirms the compact nature of the source. The absence of counterparts brighter than 27.3 mag (5σ) sets an extremely high X-ray-to-optical flux ratio (> 5260), effectively ruling out more common soft X-ray emitters such as AGN, cataclysmic variables, or the active coronae of late-type stars. At the same time, this limit is still small compared to the X-ray-to-optical flux ratios of XINSs (e.g. the optical counterparts of the predominantly thermally emitting ‘magnificent seven’ imply log(f_X/f_opt)≥4; see Fig. 7 in Kurpas et al. 2024a). Deeper optical observations of J0657 will be necessary to uncover its optical counterpart.

At X-ray energies, the XMM-Newton and NICER observations have both unveiled the spin period of the neutron star at approximately 261 ms (with a pulsed fraction of (13.7 ± 2.8)% in the EPIC pn data). The coherent timing analysis of the two epochs defines the pulsar spin-down rate at $\dot{P}=6_{-4}^{+11}\times {10}^{-15}$ s s⁻¹. However, further observations and a longer time baseline are needed to unambiguously determine a timing solution.

Fig. 6. Diagram of spin-down rate vs. spin period of the pulsars from the ATNF catalogue (version: 2.1.1; Manchester et al. 2005). The location of J0657 is indicated as a blue data point with error bars. We highlight the position of the ‘three musketeers’ (PSR B0656+14, Geminga, and PSR B1055–52; magenta), high magnetic field spin-powered pulsars (red; same selection as in Kurpas et al. 2024a), the ‘magnificent seven’ (green), and magnetars (orange).

In studying the phase-integrated X-ray spectrum of the source, we found that a PL continuum with a steep spectral slope of 4-5 or thermal models (two BB components at kT₁ ∼ 90 eV and kT₂ ∼ 220 eV, or an NSA at log(T/K)∼5.8 and a hot BB at kT ∼ 250 eV) fitted the continuum of J0657 well. In all cases, an absorption feature around ∼0.3 keV with a width of ∼50 eV is necessary for a good fit. This feature was detected with high significance, with false-negative and false-positive rates of 6.2% and 2.7%, respectively. However, the line correlates strongly with the column density, allowing only upper limits for N_H to be derived. Given the steep power-law slope, thermal models appear more physically plausible.

In Fig. 6 we show the source in the (P,Ṗ) diagram of INSs from the ATNF pulsar catalogue (version: 2.1.1; Manchester et al. 2005). Under the usual assumption of magnetic dipole braking in vacuum (Ostriker & Gunn 1969), the magnetic field and characteristic age of spin-down, $B_{\mathrm{dip}}=3.2\times {10}^{19}\sqrt{P\dot{P}}$ G and τ_ch = P(2Ṗ)⁻¹ s, are estimated as $1.2_{-0.6}^{+0.9}\times {10}^{12}$ G and $0.7_{-0.5}^{+1.7}$ Myr, respectively, adopting the sample point distributions in Sect. 3.1. These values are typical of a middle-aged RPP.

Comparing spin-down luminosity to X-ray luminosity can provide insights into the source’s nature. This is because highly magnetised INSs can exhibit X-ray luminosities that exceed the energy available from spin-down, unlike spin-powered pulsars. Using the magnetic dipole braking model, we estimate the spin-down luminosity of J0657 to be Ė = 6.3 × 10⁴⁶(ṖP⁻³) g cm² = $2.0_{-1.5}^{+4}\times {10}^{34}$ erg s⁻¹. The X-ray luminosity was calculated using the best-fit double blackbody model III in Table 4 and the Stefan-Boltzmann law. To estimate the distance to the source, we applied the upper limit on the hydrogen column density (N_H) and used the 3D-N_H Tool (Doroshenko 2024), resulting in a distance upper limit of < 1.5 kpc. In Fig. 7, we present the 1σ confidence regions and the positions of various thermally emitting INSs from Potekhin et al. (2020). The diagram illustrates that, even at a distance of 1.5 kpc, the absolute X-ray luminosity of J0657 is roughly 100 times less than its spin-down luminosity, placing it comfortably within the region occupied by RPPs.

Fig. 7. X-ray vs. spin-down luminosity diagram showing XINSs and pulsars with notable thermal emission (Potekhin et al. 2020). The black point with error bars represents the 1σ confidence region for J0657 at a distance of 1 kpc, based on XMM-Newton and NICER timing data (Sect. 3.1). The shaded grey area shows the 1σ region for distances between 0.1 kpc and 3 kpc. The dashed line marks the identity line.

Young RPPs, such as the Crab pulsar, predominantly emit non-thermal X-rays that overshadow thermal radiation from their cooling surfaces (see e.g. Becker et al. 2009, for an overview). Older pulsars (ages greater than a few million years) show a weak thermal component, likely from small heated polar caps, alongside the dominant power-law emission. For RPPs aged around 10⁵ − 10⁶ years, thermal emission is expected in soft X-rays (less than 1 − 2 keV), with magnetospheric emission dominating above 2 keV. The spectrum of J0657 suggests it may be an intermediately aged RPP, as it requires both cold and warm blackbody components to fit the data. This is consistent with models for the ‘three musketeers’, where, additionally, power-law tails with spectral indices α = 1.7 − 2.1 were detected at low flux levels (0.3 − 1.7% of the source luminosity; De Luca et al. 2005). Similarly, the upper limit for J0657 is (2.1 ± 1.8)%, indicating that any potential power-law tail is consistent with these faint levels.

Magnetospheric emission from RPPs may be detectable in both radio and gamma-ray wavelengths. Therefore, even if current X-ray observations are too shallow to identify non-thermal emission processes, detecting the source at other wavelengths might provide indirect evidence of such processes. Notably, the FAST observation of J0657 failed to detect any modulation, resulting in a stringent 10σ upper limit of 1.4 μJy at 1.5 GHz. This result is reminiscent of the behaviour observed in Geminga, which, despite numerous radio campaigns, has only shown sporadic claims of detection at low frequencies (< 100 MHz; see Malov et al. 2015; Maan 2015, and references therein). In contrast, Geminga stands out as one of the most prominent gamma-ray sources, being among the first ever identified (Bignami & Caraveo 1996).

Given the similarities between the X-ray properties of J0657 and Geminga, we examined Fermi-LAT data for a potential gamma-ray counterpart but found no significant spatial or pulsed signals over sixteen years of observations. We note that the X-ray timing solution from Sect. 3.1 was extrapolated over the entire time baseline of the gamma-ray analysis, as there were few, if any, gamma-ray photons that could plausibly be attributed to J0657. At the most likely distance range of 0.7 − 1.5 kpc (Table 4), the 95% c.l. upper limit on the photon flux corresponds to a gamma-ray luminosity of L_γ < (0.4 − 1.8)×10³² erg s⁻¹ in the 0.1 − 100 GeV energy range, assuming isotropic emission. This range lies at the lower end of the distribution for non-recycled pulsars with similar spin-down power (Smith et al. 2023).

The corresponding luminosity upper limit is illustrated in Fig. 8, alongside the 274 catalogued pulsars, including millisecond (P < 30 ms), radio-loud (S₁₄₀₀ > 30 mJy, where S₁₄₀₀ is the radio flux density at 1400 MHz) and radio-quiet gamma-ray pulsars and their upper limits. The XINS’s efficiency in converting spin-down power to gamma-ray luminosity, η ≡ L_γ/Ė < 0.9%, is also notably lower than most of non-recycled pulsars with similar Ė. The position of the high-latitude X-ray pulsar Calvera (1RXS J141256.0+792204; e.g. Halpern et al. 2013; Mereghetti et al. 2021, and references therein) is also indicated. Estimated to be at a distance of 3.3 kpc, Calvera is a spin-powered pulsar that may have been born at a significant height above the Galactic disk. Despite its relatively high Ė, it does not exhibit radio or gamma-ray emissions. Given these characteristics, it is speculated that Calvera could be a descendant of the central compact object class – a group of enigmatic, radio-quiet young neutron stars found in supernova remnants (see e.g. De Luca 2017).

Fig. 8. Gamma-ray luminosity (0.1 − 100 GeV) vs. spin-down power of pulsars from the Third Fermi-LAT catalogue (adapted from Smith et al. 2023). The black diamond with error bars and the shaded grey area represent the 95% c.l. upper limit on the energy flux of the target, 6.8 × 10−13 erg s−1 cm−2, at a distance range of 0.7 − 1.5 kpc. The positions of the ‘three musketeers’ and Calvera are highlighted (see legend). The upper diagonal line shows where 100% of spin-down power is converted to gamma-ray flux, while the lower diagonal line represents the heuristic relation $L_{\gamma }\propto \sqrt{\dot{E}}$, for reference.

A potential explanation for the absence of magnetospheric emission in certain INSs is the misalignment of the beamed non-thermal emission with the observer’s line of sight. This idea has been proposed to account for the lack of detectable radio emission in nearby, long-spin-period XINSs (Kondratiev et al. 2009). Nearly half of non-recycled gamma-ray pulsars are radio-faint (70 out of 150), and 95% of radio pulsars with similar Ė remain unseen at GeV energies (see Smith et al. 2023, for details). While the narrow radio beam is typically aligned with the INS polar regions, the difference between radio and gamma-ray emissions suggests that gamma-rays originate from a different region in the outer magnetosphere. Models (e.g. Cheng & Zhang 1998; Takata et al. 2006, 2008) propose that gamma-rays cover a wider angle than the radio beam, which could explain the phenomenon of radio-quiet pulsars among RPPs.

In addition, both beams may miss the observer. The population synthesis simulations of Johnston et al. (2020), although applied to the young and energetic portion of the spin-powered pulsar population with Ė > 10³⁵ erg s⁻¹, suggest that such a population of ‘invisible’ pulsars could account for 25% of the total. Alternatively, the absence of non-thermal emission may indicate unusual properties in the neutron star’s magnetosphere, rather than merely an unfavourable viewing geometry. Romani et al. (2011) proposed a class of sub-luminous gamma-ray pulsars with Ė > 10³⁴ erg s⁻¹ and d < 2 kpc, characterised by low gamma-ray efficiency, possibly due to a lower-altitude emission component. The XINS J0657, the ‘musketeer’ PSR B0656+14, and Calvera all belong to this category.

Absorption features in various INSs are typically linked to cyclotron resonances of charged particles in the neutron star’s magnetic field (Staubert et al. 2019), bound-bound or bound-free transitions in the atmosphere (e.g. van Kerkwijk & Kaplan 2007), or inaccuracies in temperature modelling (Viganò et al. 2014). If cyclotron resonances are the cause, the magnetic field strength for the first harmonic is given by $B_{\mathrm{cyc}}=(z+1)\frac{\mathit{mE}}{ħq}$, where m is the particle mass, E is the energy, and q is the charge. Using this formula with the absorption feature’s best-fit energy, we find $B_{\mathrm{cyc},{\mathrm{p}}^{+}}=6.5_{-0.6}^{+1.1}\times {10}^{13}$ G for protons and $B_{\mathrm{cyc},{\mathrm{e}}^{-}}=3.{56}_{-0.28}^{+0.6}\times {10}^{10}$ G for electrons. These values differ from the typical 10¹² G seen in intermediately aged RPPs, as indicated by spin-parameter limits. The higher field strength for proton cyclotron resonances indicates a more magnetised region, likely near the INS surface. In contrast, the lower field strength for electron cyclotron resonances suggests a region farther from the surface, where the magnetic field is weaker, as has been proposed for PSR B0656+14 (Arumugasamy et al. 2018).

Alternatively, the absorption feature might be due to bound-bound or bound-free transitions in the INS atmosphere or surface. If we assume that the line is caused by hydrogen ionisation, the derived magnetic field strength is $3.1_{-0.9}^{+2.8}\times {10}^{13}$ G, which exceeds both timing constraints and typical values for RPPs. However, transitions in elements with higher atomic numbers could explain the observed line, as these can occur at higher energies than hydrogen (Mori & Ho 2007; Medin et al. 2008). We also note that deeper X-ray observations might reveal additional features if transitions in higher Z elements are present.

Finally, the absorption feature might result from inadequate modelling of the surface temperature distribution. Viganò et al. (2014) suggested that this could explain similar features in some INSs, typically seen between two adjacent blackbody components with temperatures differing at least by a factor of 1.5 to 2. While this applies to the BB components in fit III (Table 4), the observed line energy is too low to fall between these components (Wien’s law suggests the BB peaks at ∼0.5 keV and ∼1.1 keV, respectively). This discrepancy might indicate the presence of additional emission components at very soft X-rays or far-ultraviolet, which would require further observations.

Additional insights into the nature of the line could be gained from phase-dependent studies. Cyclotron resonances, for instance, may show phase-dependence, as suggested for a feature in PSR B0656+14 (Arumugasamy et al. 2018; Schwope et al. 2022). However, the EPIC pn observation only permits splitting the rotation into two phase bins, and the contamination in the NICER data is too strong for a reliable phase-resolved analysis. Statistically, phase-resolved spectra fit well without an absorption component (see Table 5), and including the component yields similar line parameters to those in the phase-averaged case. Due to the limited photon count, uncertainties in the line parameters are too large to detect phase-evolution. Therefore, higher signal-to-noise X-ray data with better timing resolution are needed to investigate J0657’s phase-resolved properties.

5. Summary and conclusions

The emission properties of eRASSU J065715.3+260428 (PSR J0657+2604) are consistent with those of a middle-aged spin-powered pulsar. However, it is rare within the observed population to find pulsars of this type without detectable pulsed emission in both radio and gamma-ray energies. While this absence is most likely due to an unfavourable viewing geometry, the source’s rarity makes it particularly valuable for studying magnetospheric processes and understanding potential deviations in emission geometry. Observations at low radio frequencies could help determine whether mechanisms similar to those in Geminga and other radio-quiet gamma-ray pulsars are at play. Further X-ray observations could also provide key insights: by enabling searches for faint non-thermal emission, refining the timing solution to constrain the braking rate (critical for gamma-ray pulsation searches), and facilitating a more reliable phase-resolved spectral analysis to probe the origin of the absorption feature. Ultimately, the discovery and detailed study of this XINS highlight the unique potential of X-ray observations to uncover rare sources that conventional pulsar surveys often miss.

Acknowledgments

The authors would like to express their thanks to Wu Xiang Ping and Peng Jiang for the prompt scheduling of the FAST DDT observation. We thank Sergei Popov for valuable discussions regarding the nature of our target during the 2024 XMM-Newton Workshop. We also thank the referee, David Smith, for his suggestions for improvement, and the Fermi Helpdesk, particularly Nestor Mirabal, for their extensive guidance while analysing LAT observations. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 414059771. AMP acknowledges the Innovation and Development Fund of Science and Technology of the Institute of Geochemistry, Chinese Academy of Sciences, the National Key Research and Development Program of China (Grant No. 2022YFF0503100), the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDB 41000000), and the Key Research Program of the Chinese Academy of Sciences (Grant NO. KGFZD-145-23-15). IT gratefully acknowledges the support by Deutsches Zentrum für Luft- und Raumfahrt (DLR) through grant 50 OX 2301. This research has made use of data and/or software provided by the High Energy Astrophysics Science Archive Research Center (HEASARC), which is a service of the Astrophysics Science Division at NASA/GSFC. This work has used the data from the Five-hundred-meter Aperture Spherical radio Telescope (FAST). FAST is a Chinese national mega-science facility, operated by the National Astronomical Observatories of Chinese Academy of Sciences (NAOC). Based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme ID 111.259R.001 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. This work made use of Astropy (http://www.astropy.org): a community-developed core Python package and an ecosystem of tools and resources for astronomy (Astropy Collaboration 2013, 2018, 2022). This research made use of ccdproc, an Astropy package for image reduction (Craig et al. 2017). We derive posterior probability distributions and the Bayesian evidence with the nested sampling Monte Carlo algorithm MLFriends (Buchner et al. 2014; Buchner 2019) using the UltraNest (https://johannesbuchner.github.io/UltraNest/) package (Buchner 2021).

References

All Tables

All Figures

	Fig. 1. Lomb-Scargle periodogram for the NICER observation. The horizontal green dashed lines show the power corresponding to significance levels of 1σ to 8σ.
In the text

	Fig. 2. X-ray pulse profiles from NICER (top) and XMM-Newton (top), folded at the 261 ms modulation (see Table 3 and the text, for details). The two phase bins that were defined for the phase-resolved analysis of the XMM-Newton data are marked by the light- and dark-grey regions in the lower panel.
In the text

	Fig. 3. Corner plot depicting the parameter distribution of the coherent timing search. For presentation purposes, we subtracted the period P0 = 261.085 ms from the period values of the distribution. Vertical lines mark the 15.8%, median and 84.1% percentiles.
In the text

	Fig. 4. EPIC pn spectrum of J0657 along with the best-fit single blackbody (red, I), two blackbody (blue, II) or two blackbody with a line (green, III) models. See Table 4 for the fit parameters.
In the text

	Fig. 5. FORS2 R-band image of the field of J0657. The eROSITA (Kurpas et al. 2023, 2024a) and XMM-Newton X-ray sky positions are marked by the blue and green circles, respectively (radius marks the 1σ positional accuracy). Brown circles (with arbitrary radius of 1″) mark field sources identified from an SExtractor run.
In the text

Fig. 6. Diagram of spin-down rate vs. spin period of the pulsars from the ATNF catalogue (version: 2.1.1; Manchester et al. 2005). The location of J0657 is indicated as a blue data point with error bars. We highlight the position of the ‘three musketeers’ (PSR B0656+14, Geminga, and PSR B1055–52; magenta), high magnetic field spin-powered pulsars (red; same selection as in Kurpas et al. 2024a), the ‘magnificent seven’ (green), and magnetars (orange).

In the text

Fig. 7. X-ray vs. spin-down luminosity diagram showing XINSs and pulsars with notable thermal emission (Potekhin et al. 2020). The black point with error bars represents the 1σ confidence region for J0657 at a distance of 1 kpc, based on XMM-Newton and NICER timing data (Sect. 3.1). The shaded grey area shows the 1σ region for distances between 0.1 kpc and 3 kpc. The dashed line marks the identity line.

In the text

Fig. 8. Gamma-ray luminosity (0.1 − 100 GeV) vs. spin-down power of pulsars from the Third Fermi-LAT catalogue (adapted from Smith et al. 2023). The black diamond with error bars and the shaded grey area represent the 95% c.l. upper limit on the energy flux of the target, 6.8 × 10−13 erg s−1 cm−2, at a distance range of 0.7 − 1.5 kpc. The positions of the ‘three musketeers’ and Calvera are highlighted (see legend). The upper diagonal line shows where 100% of spin-down power is converted to gamma-ray flux, while the lower diagonal line represents the heuristic relation $L_{\gamma }\propto \sqrt{\dot{E}}$, for reference.

In the text