Detection of pulsed X-ray emission from the isolated neutron star candidate eRASSU J131716.9-402647

The X-ray source eRASSU J131716.9-402647 was recently identified from observations with Spectrum Roentgen Gamma (SRG)/eROSITA as a promising X-ray dim isolated neutron star (XDINS) candidate on the premise of a soft energy distribution, absence of catalogued counterparts, and a high X-ray-to-optical flux ratio. Here, we report the results of a multi-wavelength observational campaign with XMM-Newton, NICER and the FORS2 instrument at the ESO-VLT. We found in both the XMM-Newton and NICER data that the X-ray emission is strongly pulsed at a period of $12.757$ s (pulsed fraction $p_\mathrm{f} = (29.1 \pm 2.6)$% in the 0.2-2 keV band). The pulse profile is double-humped, and the pulsed fraction increases with energy. The XMM-Newton and NICER epochs allow us to derive a 3$\sigma$ upper limit of $\dot{P}\leq 8\times 10^{-11}$ s s$^{-1}$ on the spin-down rate of the neutron star. The source spectrum is well described by a purely thermal continuum, either a blackbody with $kT\sim95$ eV or a magnetised neutron star atmosphere model with $kT \sim 35$ eV. Similarly to other thermally emitting isolated neutron stars, we found in either case strong deviations from the continuum, a broad absorption feature at energy $\sim260$ eV and a narrow one around $590$ eV. The FORS2 instrument at ESO-VLT has not detected the optical counterpart ($m_\mathrm{R}>27.5$ mag, $5\sigma$ detection limit), implying an X-ray-to-optical flux ratio of $10^4$ at least. The properties of eRASSU J131716.9-402647 strongly resemble those of a highly magnetised isolated neutron star and favour an XDINS or high-B pulsar nature.


Introduction
The seven known X-ray dim isolated neutron stars (XDINSs) are famous for their thermal spectra, which are seemingly devoid of accretion or magnetospheric activity, have long spin periods ranging from 3 s to 17 s, and have higher magnetic field strengths (10 13 -10 14 G) than rotation-powered pulsars of similar ages (Turolla 2009).Interestingly, a more complex surface temperature distribution and high energy excess has been detected in their cumulative X-ray spectra only recently (Yoneyama et al. 2017(Yoneyama et al. , 2019;;Dessert et al. 2020;De Grandis et al. 2022).Although the currently observed sample is small, these objects might constitute a substantial part of the Galactic population of isolated neutron stars (INS).This is implied from their estimated birthrate, which could be on par with that of the much more numerous class of rotation-powered pulsars (Keane & Kramer 2008).
The seven confirmed XDINSs have all been identified in Roentgen Satellit (ROSAT) observations (Haberl 2007).Additional searches with ROSAT and the Neil Gehrels Swift Observatory, the ever-growing dataset from the pointed X-ray Multi-Mirror Mission (XMM-Newton), and the new all-sky survey at X-ray energies conducted with the extended Roentgen Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and NASA (observation 0921280101).
One of the recently identified candidates is the X-ray source eRASSU J131716.9-402647(hereafter dubbed J1317; Kurpas et al. 2023).Detected in the eROSITA all-sky surveys, with a predominantly thermal X-ray spectrum that is best fit by a single absorbed blackbody model with an effective temperature of 110 eV and a high X-ray-to-optical flux ratio that implies a compact nature, it was found to fit well into the known XDINS population.However, the data available at that time did not allow ruling out the existence of non-thermal emission.Similarly, the low photon count and survey exposure gaps 1 prevent a meaningful timing analysis.Therefore, an unambiguous classification of the source has not been possible yet, as is the case for other proposed XDINS candidates.To shed light on its nature, we carried out follow-up observations of J1317 with the Neutron star Interior Composition Explorer (NICER; Gendreau et al. 2016), XMM-Newton (Jansen et al. 2001), and the Focal Reducer/low dispersion Spectrograph 2 (FORS2) instrument at the ESO Very Large Telescope (ESO-VLT; Appenzeller et al. 1998).Notes.For XMM-Newton, we provide the EPIC-pn values in the table. (a) Modified Julian date relative to the start of the observation. (b) Where applicable, the exposure times are free from periods of high background activity, but they are uncorrected for vignetting and CCD dead time. (c) Total X-ray photon events detected in the 0.3-2.0keV energy band.
In this paper, we report the results of this observational campaign and discuss the implications for the nature of the source.We start in Sect. 2 by describing the new data sets and reduction steps.The results of the X-ray spectral and timing analysis and optical imaging are presented in Sect.3. Our conclusion and outlook are outlined in Sect. 4.

X-ray observations
The X-ray source J1317 was observed by NICER for 49 ks in March and April 2023 (Table 1).The raw data files were reduced with the NICERDAS software tools (Version: 2023-08-22_V011a) that are distributed with the HEASoft release 6.32.1.We extracted the event lists using the nicerl2 task and adopted a conservative event screening by setting the "saafilt" variable to "yes" and "saafilt_nicer" to "no", the "underonly_range" was set from 0-50, the "overonly_range" from 0-5, "cor_range" to "1.5-" and "elv" to a value of 30.Verifying that no periods of high background remained in the single observation files, we merged all observations with niobsmerge.We extracted the spectrum with the nicerl3-spect pipeline and used the SCORPEON2 model for the background estimation.Similarly, the nicerl3-lc pipeline was used to extract light curves.We applied the barycentric correction to the 55 375 cleaned and filtered events (0.3-2 keV band) with the barycorr task (ephemeris: JPLEPH.405,using the improved XMM-Newton based sky position).
The XMM-Newton observation of J1317 was performed in July 2023 for 36 ks (Table 1).The European Photon Imaging Camera's (EPIC) pn and MOS detectors, equipped with the THIN filter, were operated in full frame and small window mode (73.4 ms and 300 ms time resolution3 , respectively).We extracted data products using the Science Analysis System (SAS; Version: 21.0.0) and found the observation to be free from periods of enhanced background activity.To determine the source position, we performed source detection with the task edetect_stack in the standard XMM-Newton energy bands 0.5-1 keV, 1-2 keV and 2-4.5 keV, using the exposures from both the EPIC-MOS and pn detectors.We then applied the task eposcorr to refine the astrometry.We compared the position of all detected X-ray sources with a detection likelihood value "EP_DET_ML" above 150 with those of objects in the Guide Star Catalog (Version, 2.4.2;Lasker et al. 2008), whose positional errors are below 1.5 .The number of considered matches, the applied offset in RA and Dec from the initial XMM-Newton source detection localisation, and the updated source position are shown in Table 2.The results of the source detection were used to define a background region on the same CCD chip as the target devoid of other X-ray sources.The optimal centring and extraction region for the target was determined with the task eregionanalyse in the 0.2-2.0keV band.We decided to truncate the optimal size of 40 to a value of 16 to minimise contamination by a neighbouring hard X-ray source located 26 north of the target (Fig. 1).The event lists were further cleaned by applying the recommended filters "FLAG = 0", "PATTERN ≤ 4" for EPIC-pn and "PATTERN ≤ 12" for EPIC-MOS.Likewise, we followed the XMM-Newton data analysis threads4 to extract the spectra.For the timing analysis, the 10 374 source events from EPIC-pn (0.2-2.0 keV) were barycentric corrected with the barycen task (ephemeris: DE-405) using the updated source position.

Optical observations
Follow-up ESO-VLT observations of the field of J1317 were carried out on May 12 and 18, 2023 (Table 1).We used the FORS2 instrument equipped with the R_SPECIAL filter to obtain deep imaging under generally photometric conditions.Three exposures affected by worsening seeing conditions (>1.3 ) on May 12 were not included in the analysis.We applied the standard FORS2 reduction pipeline within the ESO Reflex environment (Freudling et al. 2013) to correct the images for bias and flat field.The individual science exposures were astrometrically corrected with astrometry.net(Lang et al. 2010).To stack the images, they were aligned with the task wcsalign (Astropy Collaboration 2013Collaboration , 2018Collaboration , 2022)).

Timing analysis
In comparison to the science modes of XMM-Newton we used, the time resolution of NICER allows us to search a broader frequency range for periodic modulations.We  (Lomb 1976;Scargle 1982) for spin periods between 0.01 s-60 s (Fig. 2).The key parameters of the search we conducted are summarised in Table 3.The periodogram shows a significant periodic signal and its first harmonic at a period of 12.76 s and 6.37 s, respectively.The existence of this modulation is confirmed by the XMM-Newton data because the Lomb-Scargle periodogram computed in the range of 0.3-60 s, using photons from EPIC-pn binned into a 150 ms light curve, also detects the modulation with high significance.
XMM-Newton allows us to verify that the nearby but fainter Xray emitting source north of the target position (see image in Fig. 1) is not the origin of the pulsed emission.A corresponding search over the XMM-Newton light curve of the source did not reveal any significant periodic signal because the highest peak possesses a false-alarm probability very close to 1.
To improve the period estimation and derive its significance, we applied the Bayesian-based folding method of Gregory & Loredo (1996).The method calculates a prescription of frequency-dependent odds ratio O m that the data favour a periodic model with m phase bins over the null (unpulsed) model.Following Zavlin et al. (2000), we adopted a characteristic number of m max = 12 and a frequency interval of 20 µHz centred on the fundamental frequency ν 0 = 0.078388 Hz.The resulting odds ratio values of O per = 9 × 10 29 and 5 × 10 48 , obtained for NICER and XMM-Newton, respectively, give a probability of virtually 100% that a periodic signal is present in the tested frequency range (see Table 3 for details).The corresponding 68% confidence intervals are P NICER = 12.757129(16) s and P XMM = 12.75707(7) s.The accuracy on the period is not high enough to phase-connect the two data sets because the more accurate NICER timing solution predicts 683 991.2 ± 0.9 spin cycles between the time of phase zero of both observations.Nonetheless, the current estimates rule out spin-down values in excess of Ṗ ≤ 8 × 10 −11 s s −1 , at the 3σ confidence level.
With the corresponding best period estimates, we folded both light curves in phase and chose the minimum of the pulse profile as the time of phase zero (see Table 3 and Fig. 3).Consistent with the results of the Lomb-Scargle the pulse profile is clearly double-humped.Interestingly, the minima between the pulses are of different depth and the pulses are of different height, the second being weaker than the first.We derived pulsed fractions of (10.7 ± 1.2)% (0.3-2.0 keV) for NICER and (29.1 ± 2.6)% (0.2-2.0 keV) for XMM-Newton, respectively.The background was not subtracted for either value.Thus, the discrepancy may be explained by a stronger background event contamination for NICER.In Fig. 4, we show the energy and phase dependence of the pulse for EPIC-pn.The pulse at phase 0.7 seems to be strongest at intermediate to high energies, whereas the pulse around phase 0.25 seems to be strongest at soft energies.Similarly, the pulsed fraction increases with energy from 15(7)% (0.2-0.3 keV) to 65( 14)% (1-1.5 keV).

X-ray spectral analysis
We used XSPEC (version 12.13.0,Arnaud 1996) to fit spectral models to the XMM-Newton and NICER data sets.The spectra were binned with at least 25 counts per spectral bin to allow for the use of the χ 2 statistic.We accounted for the interstellar absorption by applying the tbabs model with elemental abundances described in Wilms et al. (2000).When spectra from several instruments (EPIC-pn, EPIC-MOS1, EPIC-MOS2, and NICER) were fitted simultaneously, we added a constant factor to account for calibration uncertainties.When fitting the NICER spectrum, the background was simultaneously fitted adopting the SCORPEON model.We found that all spectra are affected by the hard source located north of the target.To account for the level of contamination and estimate its dependence with energy, we extracted the spectrum of the hard source from the EPIC-pn exposure.We used a small aperture of 6 to prevent bias from the brighter point spread function (PSF) of the target.Due to the low count number (357 events extracted in the 0.2-12 keV energy range), we adopted the Bayesian X-ray Analysis (BXA; A164, page 3 of 10 Table 3. Timing parameters from the NICER and XMM-Newton observation.

NICER
XMM-Newton N photon 55 375 10 374 Energy range [keV] 0.3-2.00.2-2.0∆T (a)  [s] 246 057 36 349 ν LS,min -ν LS,max [Hz] 0.0167-100 0.0167-3.33LS significance (b)  5.4σ 6.7σ 0.07838754(10) 0.0783879(5) Period [s] 12.757129( 16) 12.75707(7) Pulsed-fraction (d) [%] 10.7 ± 1.2 29.1 ± 2.6 Time-of-phase zero (e) 60 032.590341 60 133.582978 Notes. (a) Time between first and last photon. (b) For the fundamental peak, we show the significance that the observed Lomb-Scargle power is not due to white noise.The chance that the observed peak is due to noise was computed from the false-alarm probability (Baluev 2008). (c) Probability of detection of a periodic signal in the frequency interval we probed.Computed via O per /(1 + O per ), with Zavlin et al. 2000, for more details). (d) Computed from the maximum and minimum light-curve count-rate values (R) via Rmax+R min (Fig. 3). (e) Reference time in MJD used to align the NICER and XMM-Newton pulse profiles.Two cycles are shown for visual purposes.Buchner et al. 2014) package to fit the unbinned data.The spectrum of the contaminant is well described by an absorbed powerlaw model with a photon index of Γ = 1.81 +0.15 −0.13 , N H = (10 ± 4) × 10 20 cm −2 and f X = 8.4(5) × 10 −14 erg s −1 cm −2 (0.2-10 keV), typical of that of an active galactic nucleus (AGN).The column density is consistent with the Galactic value of 6.7 × 10 20 cm −2 in the direction of the source (HI4PI Collaboration 2016).We added this power-law model as an additional background component to fits of the NICER spectrum.For XMM-Newton the contamination can be mitigated if the energy range is constrained to 0.2-2.0keV.We found XMM-Newton and NICER to give consistent results when NICER is limited to the energy range of 0.4-10 keV.We list the fit results of all models we used in Table 4.The number of detected photons in EPIC-pn, normalised by the median value derived from the sum of all photons that were detected in a certain energy bin column, is colour-coded.At the top, we indicate the pulsed fraction, determined from events in a certain energy bin, whereas to the right, we plot the normalised light curves in the energy bands 0.2-0.4keV (blue), 0.4-0.8keV (orange), and 0.8-2 keV (green).They were computed by summing up all events in a certain phase bin and energy range and normalising again by the median number of counts in the given energy range and all phase bins.
We began the spectral analysis by fitting a single blackbody (BB), neutron star atmosphere (NSA; Pavlov et al. 1995;Zavlin et al. 1996) or power-law (PL) model to the EPIC-pn spectrum.This approach is motivated by the fact that the spectral analysis of the eROSITA data could not exclude any of these models (Kurpas et al. 2023).For all three models, we found high systematic residuals, a poor fit quality, and column densities in excess of the Galactic value (see Table 4 and Fig. 5, for the BB case).The addition of a PL component on the thermal models (BB, NSA), as usually observed as a hard tail in the spectrum of middle-aged spin-powered pulsars (e.g.Schwope et al. 2022), does not improve the fit.Combining multiple BB components, which is indicative of a more complex surface temperature distribution, did not improve the fit either.Interestingly, we found that the addition of a multiplicative broad Gaussian absorption line (GABS model) at the low-energy end of the spectrum (200−400 eV) significantly improved the fit quality.For a PL continuum (combined with a single GABS component), the resulting fit shows a good agreement between data and model and no significant residuals that might imply the necessity for additional model components.The fit of the thermal models (both BB and NSA) with one GABS component, although an improvement with respect to the single-component case, still shows high χ 2 ν values and systematic residuals consistently around energies of 550−600 eV (see Fig. 5, again for the BB case).The addition of a second GABS component improves the results in both cases, but it introduces degeneracy between N H and the parameters of the lines.Fixing the widths of the two Gaussian components to their best-fit values, σ 1 = 150 eV and σ 2 = 15 eV, leads to slightly more accurate line energies and smaller errors on the N H (shown for the BB case in Table 4), but does not allow us to constrain the exact value of the column density.Thus, the N H values we obtained can only be regarded as an upper limit, but they strongly favour a Galactic source.
For a BB continuum, we studied the significance of the detected lines by computing the false-negative rates (an existing feature is not identified) and false-positive rates (the fit A164, page 4 of 10  PL  5)  Notes.We give the 1σ confidence intervals for the estimated parameters. (a) We assumed a 1 kpc distance for the blackbody emission radius at infinity. (b) The equivalent width (EW) was estimated via Fc−Fo Fc dE, with F c being the continuum and F o the observed flux.The errors state the maximum and minimum EW values obtained from all possible combinations of the upper and lower 1σ confidence interval limits of the model parameters. (c) The absorbed model flux covers the 0.2-10 keV range. (d) The model parameters assume a canonical neutron star with 1.4 M and 12 km radius.statistics are improved only by chance when absorption components are added).We used the best-fit BB2GABS model and the XSPEC fakeit command to simulate 1000 spectra and fitted them with a single BB, BBGABS, and BB2GABS model.We found that a BB model was never preferred over the BBGABS fit and that a BBGABS only gave better fits than a BB2GABS model in 2% of the cases.Since the line parameters were generally well recovered, the simulations imply very low false-negative rates close to 0% for the broad feature at ∼260 eV and 2% for the narrow line at ∼590 eV.For the false-positive case, we simulated 1000 BB and BBGABS spectra each and then conducted fits with the simulated model and with one model that contained an additional absorption component.In both cases, we found that similarly large improvements in the fit statistic (∆χ 2 ν = 2.4 for the broad feature and ∆χ 2 ν = 0.65 for the narrow one) by the addition of a line were never observed in the simulations.This implies falsepositive rates close to 0% for both features.We conclude that both features seem to be detected with very high significance.
We next tried fits with other absorption feature components available in XSPEC (GAUSS and edge), combining them with a BB continuum.Fits with the GAUSS component reproduce the results from the GABS models.On the other hand, three edge components are necessary to properly fit the spectrum (χ 2 ν (ν) = 1.18(20)), suggesting a possible third line at 864 +22 −25 eV.However, we note that the resulting N H = 14 +4 −3 × 10 20 cm −2 significantly exceeds the Galactic value and that the fit statistics is worse than in the GABS and GAUSS cases.
The NSA model fits have shown that in comparison to nonmagnetised models, better χ 2 ν values are obtained when the magnetic field is taken into account.They allow us to obtain distance estimates.We do not list them in Table 4, but found the largest distances for magnetised NSA models with two absorption lines, converging to 100 +70 −16 pc (10 12 G) and 130 +160 −18 pc (10 13 G), respectively.These may be comparable (within 2-3σ) to those inferred for the known population of XDINSs (see e.g.Motch et al. 2009, and references therein), but could be too small considering the much lower flux of J1317.
Alternatively, composite models consisting of two blackbody components or of a blackbody plus power-law, in both cases modified by a Gaussian absorption feature (2BBGABS and BBPLGABS, respectively), also fit the data well (χ 2 ν (ν) = 1.02(21) A164, page 5 of 10 and χ 2 ν (ν) = 0.87( 21)).However, for both models, the column density is a factor of 1.5-5 in excess of the Galactic value (considering the 1σ confidence interval).For the 2BBGABS case, the BB radius at 1 kpc distance is too large for an INS nature.Assuming a canonical radius of 12 km, this would imply distances between 40 and 400 pc that are in stark contrast to the high N H values.In the case of the BBPLGABS model, the photon index of the power-law component is unreasonably steep (Γ = 12.4 +1.6 −2.8 ), while the thermal component only dominates the high-energy end of the spectrum.This is at odds with what is observed in, X-ray binaries and rotation-powered pulsars, for instance, making these composite models overall less compelling than those with multiple lines.
The investigation of a hard excess in the spectrum of J1317 is made difficult at the current signal-to-noise ratio by contamination from the neighbouring X-ray source.In order to constrain it, we extracted the EPIC-pn spectrum using a large extraction region radius of ≈40 , which encompasses PSF counts from both sources.We then fitted the spectrum in the energy range of 0.2-8 keV.We defined a model consisting of an absorbed BB component with two absorption lines and an independently absorbed PL component, with both the column density and photon index parameters fixed to the best-fit values obtained for the contaminant.The fit is virtually identical to the BB2GABS fit in Table 4, with no indications of a systematic excess emission at harder energies.The BB and PL components converge to flux values of 4.62(5) × 10 −13 erg s −1 cm −2 and 8.8(6) × 10 −14 erg s −1 cm −2 (both 0.2-10 keV), respectively.These values agree well with those obtained from the single fits to J1317 and the contaminant spectrum.Notes.We give the 1σ confidence intervals for the estimated parameters.The N H is set to 5×10 20 cm −2 , the line widths are set to 150 eV and 15 eV for the first and second line, respectively.The fit converged to a temperature of kT = 91.8+1.5 −1.5 eV and fit statistic of χ 2 ν (ν) = 1.01( 72). (a) We assumed a 1 kpc distance for the blackbody emission radius at infinity. (b) The equivalent width (EW) was estimated the same way as for the phase-averaged fits.See notes of Table 4 for more details.
Using the period estimate from XMM-Newton (Table 3), we performed phase-resolved spectroscopy with the EPIC-pn dataset by splitting the photons into four phase bins (0.875-0.125, 0.125-0.375,0.375-0.625,and 0.625-0.875)including 2000-3000 photons each.The four spectra were then simultaneously fitted with a BB, BBGABS, and BB2GABS model.To lift the degeneracies and aid the fits, the N H was fixed to the phaseaveraged upper limit of 5×10 20 cm −2 , the line widths of the broad and narrow feature were fixed to 150 eV and 15 eV, respectively, and the effective temperature was tied to the same value for all phases.We found that a single BB model does not give a convincing fit (χ 2 ν (ν) = 9.60( 82)), but the fit can be improved when broad lines at ∼300 eV are included at all phases (χ 2 ν (ν) = 1.38( 74)).This fit shows strong residuals at ∼590 eV that are most prevalent at phase 0.125-0.375.Including a second narrow line in this phase bin yields a near perfect fit (χ 2 ν (ν) = 1.01(72); see Table 5).Fitting this narrow feature in other phase bins gives similar line energies around 590 eV, but only marginally improves the fit statistic (χ 2 ν ∼ 0.95).We found that leaving the line widths or effective temperatures free to vary during the fits only slightly improves the fit statistic, but results in decreased precision in the parameter estimation, such that no significant changes in temperature or line width can be observed.We studied the detection significance of the features by simulating for each phase and feature 1000 spectra with and without the best-fit line.We then repeated the spectral analysis with the real and simulated data, including in each fit only one simulated spectrum at a time, to study the false-positive and false-negative rates.The line parameters of the broad features and the narrow feature at phase 0.125-0.375were generally well recovered, and all fits implied low falsepositive and false-negative rates virtually identical with zero.However, the current data do not allow us to significantly detect or exclude the existence of the narrow feature at the other phase bins because the simulations resulted in false-positive rates of 5-70% and false-negative rates of 8-47%.We conclude that the broad feature appears to be observable in all phases.The narrow feature could be phase dependent because it was only significantly detected in one phase bin.However, additional observations are necessary to fully characterise its phase-dependence.
We tried simultaneous phase-averaged fits of all the available instruments (EPIC-pn, EPIC-MOS1, EPIC-MOS2, and NICER) next.The results are generally consistent with the single EPICpn fits.The only difference is that thermal models, including only a single absorption line, already fit well (as shown for the BB case in Table 4).Systematic residuals, indicating the necessity for a second line around 600 eV, are not apparent.
A164, page 6 of 10 We attribute this to the fact that the noise added by including NICER and EPIC-MOS dominates in comparison to the residuals caused by the absorption feature.Consequently, adding a second GABS component to the spectral fit improves the fit statistic only slightly and is less well constrained than EPIC-pn.We thus accept as our final model a pure blackbody continuum of 95.1 +1.7 −2.6 eV modified by a broad line at 260 +80 −50 eV, a narrow line at 587 +8 −5 eV, and cold interstellar matter with a column density below 5 × 10 20 cm −2 (see line 3 in Table 4).

VLT observation
With respect to previous imaging with the Southern African Large Telescope (SALT; Kurpas et al. 2023), the deeper ESO-VLT observations allow us to further improve the limit on the X-ray-to-optical flux ratio for this source and to search for faint optical counterparts.We applied the SExtractor software (Bertin & Arnouts 1996) to detect all sources near J1317 and mark them in Fig. 6.The updated XMM-Newton EPIC position of the target, although consistent within 2σ, is 1.3 offset from that derived from the eROSITA data.This brings J1317 closer to two field sources that are now separated by only 1.7 and 2.5 (see sources to the east in Fig. 6).The magnitudes of these sources were measured at 26 mag and 24.5 mag, respectively.They were not detected in the SALT image presented in Kurpas et al. (2023), since their magnitudes are close to the detection limit.Nonetheless, the association of the target with any of these possible optical counterparts is unlikely, with angular separations at 7σ and 11σ (Table 2).
We used the nightly determined calibration zero points and extinction values provided by ESO 5 for the magnitude calibration of the FORS imaging.Assuming a 5σ detection limit, the equation for the magnitude limit of a point source in an optimal Gaussian aperture was used to estimate the image depth.Based on the parameters shown in Table 2, we obtained a limiting magnitude of 27.5 mag at the position of the target, resulting in an X-ray-to-optical flux ratio of 10 4 at least.Such an extreme value can only be observed among INSs (Fig. 7).

Discussion and outlook
In this work, we report the results of dedicated NICER, XMM-Newton, and ESO-VLT follow-up observations of the INS candidate J1317.This campaign unveiled a highly significant periodic signal at a period of 12.757 s, very likely associated with the spin of the neutron star.The pulse profile is clearly double-humped, and we found the pulsed fraction to increase with energy.The results are well in line with expectations for a highly magnetised INS.
In Fig. 8 we plot the spin versus spin-down diagram of neutron stars from the ATNF Pulsar Catalogue 6 (Version 1.70; Manchester et al. 2005).The detected periodicity places J1317 where we would expect to find an XDINS, magnetar or high-B pulsar (see e.g.Enoto et al. 2019, and references therein, for an overview of the Galactic INS population).The XMM-Newton and NICER epochs only allow us to derive an upper limit on the spin-down of the source, Ṗ ≤ 8 × 10 −11 s s −1 (3σ).Under the usual assumption of magnetic dipole braking in vacuum (Ostriker & Gunn 1969), the dipolar magnetic field 5 The magnitude zero points and extinction values can be accessed at: http://archive.eso.org/qc1/qc1_cgi?action=qc1_browse_ table&table=fors2_photometry 6 https://www.atnf.csiro.au/people/pulsar/psrcatFig. 7. X-ray-to-optical flux ratios of several classes of X-ray emitters as a function of hardness ratio (see Kurpas et al. 2023, for details).The updated position of the target is denoted by the lower limit in black and is compared to those of the five XDINSs on the German eROSITA footprint (western hemisphere from the Galactic centre; magenta data points).
in the equator and characteristic age of J1317 are constrained to B dip = 3.2 × 10 19 √ P Ṗ G ≤ 10 15 G and τ ch ≥ 2600 yr.A better constraint on the evolutionary state of the neutron star relies on the eventual measurement of the spin-down rate and possible detection of radio and γ-ray counterparts.
The deep ESO-VLT FORS2 observations of the field of the XDINS candidate in the R_SPECIAL band revealed two faint sources in the vicinity of J1317.The association of these two sources with the neutron star candidate is unlikely however, because the updated XMM-Newton EPIC position implies a significant separation of 7-11σ from these objects.The absence of counterparts to such an extreme X-ray-to-optical flux ratio, ≥10 4 , further supports an INS nature (Fig. 7).
A164, page 7 of 10 The spectrum of J1317 cannot be described by a singlecomponent model.Good spectral fits can be obtained if the continuum (a blackbody, power-law, or a fully ionised neutron star model atmosphere) is modified by one or two Gaussian absorption features.All these models have in common that one line is located at soft energies between 200 eV and 400 eV with a width of σ ∼ 150 eV.The second line, strictly necessary for thermal continuum fits (BB and NSA) to the EPIC-pn spectrum, is located around 590 eV with a smaller Gaussian width of σ = 10−30 eV.The lines seem to be detected with very high significance.Testing other absorption components, we found that only an edge model affects the spectral results significantly because a third feature at ∼860 eV is necessary to achieve an acceptable fit.However, we note that compared to the models using a GABS component, the fit statistic is generally worse and the N H significantly exceeds the Galactic value, making them less convincing.The temperature of the continuum spectrum is kT ∼ 95 eV in the BB case and kT ∼ 35 eV for the NSA model.We note that the low NSA model temperature may predict an optical flux that is too high, as is known from the XDINSs (Pavlov et al. 1996).Similarly, the NSA model assumes uniform temperature emission from the full neutron star surface, which might explain the low distance values from this model.Other spectral models may give statistically valid fits, but are less compelling on the basis of their parameter values (e.g. a 2BBGABS model that implies a small distance of only 40-400 pc for a typical neutron star radius of 12 km, but, despite the absorption feature, converges to column densities well in excess of the Galactic value).
The eROSITA spectra discussed in Kurpas et al. (2023) did not allow us to statistically exclude a PL nature of the source.Here, a good fit can be obtained when a PL is combined with an absorption feature, but the resulting photon index of Γ = 9.5(4) is even larger than observed for eROSITA.This value is too high because neutron star non-thermal emission components are usually observed to have Γ < 3 (Becker 2009).The column density exceeding the Galactic value may imply an extragalactic nature.The only extragalactic objects that reach very large photon indices are Seyfert I galaxies (Grünwald et al. 2023).However, the very high X-ray-to-optical flux limits and the strong detected pulsation at 12.757 s make an Seyfert I galaxy nature very unlikely.Thus, it seems that the PL can be rejected as a valid spectral model.
The obtained flux values from the EPIC-pn spectral fits and from fits using all spectra (Table 4) agree with those derived from eROSITA (Kurpas et al. 2023) and thus by extension also with ROSAT.Nevertheless, if the long-term variability of J1317 is to be probed, the influence of the contaminant needs to be considered, which was not separately detected in ROSAT, in any of the eROSITA surveys, or in the independent analysis of Kurpas et al. (2023).The best-fit contaminant power-law model indicates a flux of ∼3 × 10 −14 erg s −1 cm −2 in the 0.2-2 keV band, where the emission is most prevalent in the short eROSITA observation.This is of the same size as the 1σ eROSITA flux error.Since the eROSITA confidence interval agrees well with results from XMM-Newton and NICER, we conclude that the contaminant did not significantly affect the eROSITA flux measurement.Since the estimated eROSITA count rate is comparable with the ROSAT measurement (with a similarly short exposure), the emission of J1317 seems to have been stable over the past 30 years.
The BB fits give upper limits on the absorbing hydrogen column density and can therefore be used to compute limits on the distance to the source.We used the relations given in Güver & Özel (2009) and Cardelli et al. (1989) to compute the total extinction A(5500 Å) < 0.226 for N H < 5 × 10 20 cm −2 .We compared this limit to the total extinction in the direction of J1317 as it can be inferred from the Gaia-2MASS 3D maps of Galactic interstellar dust (Lallement et al. 2022) and found them to be almost exactly in accordance with a distance of 1 kpc.This implies a radius of 5.1 +1.2 −0.5 km and a luminosity of 2.7 +1.3 −0.6 × 10 32 erg s −1 .More constraining values can be derived from the spectral fit with fixed Gaussian line widths.Here, the limit of N H < 2 × 10 20 cm −2 implies A(5500 Å) < 0.09 and indicates a distance of only 185 pc.We then computed a BB radius of 0.96 +0.13  −0.09 km and a luminosity of 9.6 +2.5 −1.8 × 10 30 erg s −1 .We note that similarly to the NSA distance estimate, the distance of 185 pc may be too small.While the value would agree with the distance to the closest XDINS RX J1856.5-3754 (van Kerkwijk & Kaplan 2007), we note that the flux of J1317 is a factor 4-20 lower than those observed for the other XDINSs (∼2×10 −12 −1×10 −11 erg s −1 cm −2 ).Only the flux of the XDINS RX J0420.0-5022 may be comparable with J1317 (Haberl et al. 2004), but the spectrum of this source is much softer, and its distance was estimated to ∼345 pc (Posselt et al. 2007).
The fact that the thermal continuum needs to be modified by multiple absorption features to properly fit the spectrum makes J1317 very similar to XDINSs and other thermally emitting INSs.In particular, a very similar broad absorption line at 200(40) eV with a Gaussian width of σ = 139 +12 −13 eV and equivalent width (EW) of 179 +3 −59 eV has been reported in the EPICpn spectrum of RBS 1223 (Borghese et al. 2017).Interestingly, RBS 1223 also exhibits a double-humped pulse profile at a similar long spin period of 10.31 s, and it is the source with the highest pulsed fraction of the seven XDINSs (19%, 0.2-1.2keV).As observed for J1317, the pulsed fraction seems to increase with energy (Schwope et al. 2005).The other XDINSs possess features with lower EW values that agree better with the more narrow line at ∼590 eV (e.g.Turolla 2009).
Absorption features in the spectra of XDINSs are expected to arise due to proton cyclotron resonances (Staubert et al. 2019), atomic bound-bound or bound-free transitions in the neutron star atmosphere (e.g.van Kerkwijk & Kaplan 2007), the presence of highly ionised oxygen in the ISM, atmosphere, or vicinity of the neutron star (Hambaryan et al. 2009;Hohle et al. 2012), or they might be of spurious nature because multi-temperature thermal components on the surface of the neutron star might mimic such features (Viganò et al. 2014).The general consensus is that the spectra combine features with different physical origins, and thus, the two detected features in the spectrum of J1317 might generally be of different nature.When we interpret them to be caused by proton cyclotron absorption, the magnetic field of J1317 could be estimated via E cyc = 0.063 B 10 13 G × (1 + z) [keV], with z being the gravitational redshift.For a canonical neutron star with M = 1.4 M and R = 12 km, z = 0.35.Based on the line energy of 260 +80 −50 eV, a proton cyclotron feature would imply a magnetic field of 5.6 +1.8 −1.1 × 10 13 G.This would indicate that J1317 is a highly magnetised neutron star, and it agrees well with the dipolar magnetic field strengths inferred for XDINSs.The line energy of the second feature would imply a magnetic field strength of 1.255 +0.018 −0.011 × 10 14 G, indicating that it could emerge from a more magnetised structure near the neutron star surface.This configuration was proposed for narrow absorption features in other thermally emitting INSs (e.g.Tiengo et al. 2013;Borghese et al. 2017).However, these features are usually too weak to be detectable in the phase-averaged spectrum and were observed to be highly phase dependent.Given the results of the phase-resolved spectroscopy, establishing or excluding a A164, page 8 of 10 10 2 10 1 10 0 10 1 10 2 Spin period [s]  Manchester et al. 2005).The period of J1317 is marked by the vertical green line.The horizontal magenta line indicates the upper 3σ spindown limit derived from the XMM-Newton and NICER observations.phase-dependent nature of the narrow feature of J1317 is difficult, however.Alternatively, the features of J1317 might be interpreted as being produced by atomic transitions in a hydrogen atmosphere, but in this case, the feature at ∼590 eV would also need magnetic field strengths in excess of 10 14 G (see Fig. 7 in van Kerkwijk & Kaplan 2007), which is higher than the dipolar fields observed in the other XDINSs.Thus, the presence of a He or higher-Z atmosphere would be more compelling because they allow for lines at the observed energy and typical XDINS magnetic field strengths (see e.g.Pavlov & Bezchastnov 2005;Mori & Ho 2007).Absorption by the presence of highly ionised oxygen would cause features between 500−600 eV.These lines are generally narrower than the line that is observed for J1317, which might imply that it is intrinsic to the source (comparable to the interpretation of the line at ∼530 eV in the RGS spectrum of RBS 1223; Hohle et al. 2012).In the last case, the features could be artefacts caused by a non-optimal spectral modelling that does not sufficiently account for the effects of a more complicated surface temperature distribution.However, Viganò et al. (2014) reported that the spectrum of RBS 1223 is difficult to model with inhomogeneous temperature distributions.Because the spectra of J1317 and RBS 1223 are very similar, this might also hold true here.Based on the available data alone, an exact characterisation of the observed features is difficult because at the moment, multiple production mechanisms are possible for each of them.Additional observations are required to attain more stringent constraints on their parameters (e.g.phase dependence, energy, width, and strength).
The soft and seemingly purely thermal spectrum and lack of significant flux and spectral variability speak against a magnetar nature.These sources, usually discovered when undergoing bright bursts of X-ray and soft γ-ray emission, possess much harder spectra, show higher persistent X-ray fluxes and strong and complex variability.Similarly, the absence of a detectable supernova remnant at X-ray energies seems to exclude a young pulsar nature.The XDINS population is oftentimes compared to the population of high-B pulsars (see Enoto et al. 2019 2020).The luminosity limits inferred for J1317 are shown in green (for a distance of 185 pc) and black (for a distance of 1 kpc).The grey box indicates the parameter space constrained by the current data, assuming distances within 0.1-1 kpc and the spin-down limit from the timing analysis Ṗ ≤ 8 × 10 −11 s s −1 (3σ).The dashed black line indicates where both luminosities are equal in value.
rotation-powered pulsars and magnetars (we indicated high-B pulsars mentioned in Ng & Kaspi 2011;Olausen et al. 2013;and Potekhin et al. 2020 in Fig. 8), these objects possess magnetic fields of B dip = 10 13 −10 14 G.However, their emission is thought to be governed by magnetic braking and not by remnant heat and the decay of the magnetic field, as is proposed for the XDINSs.In comparison with the latter, high-B pulsars seem less efficient in converting spin-down power into X-ray luminosity, as is the case for conventional pulsars (e.g.Possenti et al. 2002).Based on the spin-down rate limit derived for J1317 (Sect.3.1), we can constrain its spin-down luminosity as Ė ∼ 6.3 × 10 46 ( ṖP −3 ) g cm 2 ≤ 2.4 × 10 33 erg s −1 at the 3σ level.For comparison, in Fig. 9 we show the X-ray luminosity of several isolated pulsars as a function of the luminosity available from spin down.The parameter space covered by our current limits on J1317 is marked by a shaded grey area, assuming that the source is at a distance of 0.1-1 kpc.Based on the current limits, both a high-B pulsar and an XDINS are possible for J1317.This makes the neutron star an interesting target for a search for its radio counterpart.On the one hand, the detection of coherent radio emission may provide an independent distance estimate from the dispersion measure and a possible constraint on the pulsar braking index via a precise timing solution.On the other, it can shed light on whether XDINSs are truly radio quiet or are nearby long-period radio pulsars for which the narrow emission beam simply does not sweep over the Earth (Kondratiev et al. 2009).In comparison to the known high-B pulsars (0.15 s ≤ P ≤ 7.74 s), the spin period of J1317 is longer.At the time of writing, J1317 is only the 13th neutron star detected with a spin period in excess of 10 s.Only the XDINS RX J0720.4-3125(Hambaryan et al. 2017) and the radio detected neutron stars GLEAM-X J162759.5-523504.3 (Hurley-Walker et al. 2022), PSR J0901-4046 (Caleb et al. 2022), PSR J0250+5854 (Tan et al. 2018), and PSR J1903+0433g (Han et al. 2021) possess longer spin periods.Combined with its bright thermal luminosity, which is already conservatively comparable to that available from spin-down, this suggests that the evolutionary path is A164, page 9 of 10 considerably affected by the super strong magnetic field of the neutron star.
To summarise, the recent observations of J1317 have clearly confirmed the INS nature of the source.The detected pulsations at 12.757 s and evidence of spectral absorption features are remarkable among X-ray pulsars and may indicate the signature of a strong magnetic field.Its overall properties are in line with an XDINS nature within the scenario that they have evolved from magnetars, but they are also similar to those of high-B pulsars for which the effects of field decay have not been as significant.The XMM-Newton and NICER observations give only a weak limit on the spin-down rate of the neutron star.A dedicated coherent timing solution is necessary to further characterise the evolutionary state of the source, that is, to determine its position in the spin versus spin-down diagram of INSs and obtain more stringent estimates on the neutron star magnetic field, characteristic age, and spin-down luminosity.Similarly, observations at radio energies might be interesting because the apparent lack of radio emission in XDINSs is one of their main defining characteristics.

Fig. 1 .
Fig. 1.RGB XMM-Newton EPIC image.Red indicates counts in the 0.2-1 keV band, green in the 1-2 keV band, and blue in the 2-12 keV band.The source and background regions used in the EPIC-pn analysis are marked by the green circle and box.The X-ray sources we used to correct the X-ray sky position of J1317 are marked by cyan circles.

Fig. 4 .
Fig.4.Dependence of the pulse profile as a function of energy and phase.The number of detected photons in EPIC-pn, normalised by the median value derived from the sum of all photons that were detected in a certain energy bin column, is colour-coded.At the top, we indicate the pulsed fraction, determined from events in a certain energy bin, whereas to the right, we plot the normalised light curves in the energy bands 0.2-0.4keV (blue), 0.4-0.8keV (orange), and 0.8-2 keV (green).They were computed by summing up all events in a certain phase bin and energy range and normalising again by the median number of counts in the given energy range and all phase bins.

Fig. 5 .
Fig. 5. EPIC-pn spectrum in the energy range 0.2-2 keV.We indicate the best model and residuals for a single absorbed BB model (red), an absorbed BB with a single absorption component (blue), and an absorbed BB with two absorption components (green).

Fig. 6 .
Fig. 6.FORS2 R_SPECIAL-band image of the field of J1317.The circles indicate the 1σ positional uncertainty estimated from XMM-Newton (green) and eROSITA (blue).We searched for possible counterparts using the SExtractor software.The identified field objects are marked in brown (circles have arbitrary radii of 1 ).

Fig. 8 .
Fig.8.Spin vs. spin-down diagram presenting the spin properties of the known isolated neutron star population.In the background, we plot all the objects collected in the ATNF pulsar database (Version 1.70;Manchester et al. 2005).The period of J1317 is marked by the vertical green line.The horizontal magenta line indicates the upper 3σ spindown limit derived from the XMM-Newton and NICER observations.

Fig. 9 .
Fig.9.X-ray luminosity as a function of spin-down luminosity for the collection of thermally emitting isolated neutron stars fromPotekhin et al. (2020).The luminosity limits inferred for J1317 are shown in green (for a distance of 185 pc) and black (for a distance of 1 kpc).The grey box indicates the parameter space constrained by the current data, assuming distances within 0.1-1 kpc and the spin-down limit from the timing analysis Ṗ ≤ 8 × 10 −11 s s −1 (3σ).The dashed black line indicates where both luminosities are equal in value.

Table 1 .
All X-ray and optical observations.

Table 4 .
X-ray spectrum fitting results.
, and references therein).Located at the intersection of conventional