Thermal X-ray emission identified from the millisecond pulsar PSR J1909-3744

Pulsating thermal X-ray emission from millisecond pulsars can be used to obtain constraints on the neutron star equation of state, but to date only five such sources have been identified. Of these five millisecond pulsars, only two have well constrained neutron star masses, which improve the determination of the radius via modelling of the X-ray waveform. We aim to find other millisecond pulsars that already have well constrained mass and distance measurements that show pulsed thermal X-ray emission in order to obtain tight constraints on the neutron star equation of state. The millisecond pulsar PSR~J1909--3744 has an accurately determined mass, M = 1.54$\pm$0.03 M$_\odot$ (1 $\sigma$ error) and distance, D = 1.07$\pm$0.04 kpc. We analysed {\em XMM-Newton} data of this 2.95 ms pulsar to identify the nature of the X-ray emission. We show that the X-ray emission from PSR~J1909--3744 appears to be dominated by thermal emission from the polar cap. Only a single component model is required to fit the data. The black-body temperature of this emission is kT=0.26\ud{0.03}{0.02} keV and we find a 0.2--10 keV un-absorbed flux of 1.1 $\times$ 10$^{-14}$ erg cm$^{-2}$ s$^{-1}$ or an un-absorbed luminosity of 1.5 $\times$ 10$^{30}$ erg s$^{-1}$. Thanks to the previously determined mass and distance constraints of the neutron star PSR~J1909--3744, and its predominantly thermal emission, deep observations of this object with future X-ray facilities should provide useful constraints on the neutron star equation of state.


Introduction
Fifty years after the discovery of neutron stars, the nature of the material making up their core remains largely unknown. At densities above a few times the saturation density of nuclear matter, certain models predict the existence of exotic components such as hyperons or unconfined quarks (e.g. Lattimer & Prakash 2007). Neutron stars constitute the only medium in which nuclei exist in extremely dense but relatively cold environments. It is therefore essential to understand the equation of state of this material (meaning its pressure as a function of energy density) if we wish to have a full understanding of the baryonic matter composing the Universe.
Measuring the mass and radius of a neutron star allows us to constrain the density and pressure of the neutron star. Tight constraints should be possible with future gravitational wave observations of neutron star mergers, but to date, only weak constraints are available (e.g. De et al. 2018;Abbott et al. 2018). Using electromagnetic observations, accurate mass measurements have been made for a number of radio pulsars, for example, 1.18 +0.03 −0.02 M , for the case of one of the pulsars in PSR J1756-2251 (Faulkner et al. 2005) or 1.4414 ± 0.0002 M , for PSR B1913+16 (Weisberg & Taylor 2005). While there are some accurate mass measurements already available, the situation is different with radii that are much harder to constrain (e.g. Özel et al. 2016a;Miller & Lamb 2016). A promising way of inferring radii (and masses) is through modelling the X-ray waveform of a neutron star with a fairly low magnetic field and which shows predominantly thermal X-ray emission from its polar caps. This provides a measure of the mass and radius of the star (e.g., Pavlov & Zavlin 1997;Bogdanov et al. 2007;Leahy et al. 2011). Millisecond pulsars have magnetic fields, B 10 9 G, which are weak enough that they do not modify the opacity of the neutron star atmosphere (e.g. Heinke et al. 2006). Five millisecond pulsars have been shown to have predominantly thermal X-ray emission and show X-ray pulsations. These are PSR J0030+0451, PSR J2124-3358, PSR J0437-4715, PSR J1024-0719 and PSR 1614-2230 (e.g., Zavlin 2006;Bogdanov et al. 2008;Bogdanov 2013;Pancrazi et al. 2012). Indeed, loose constraints have been obtained for the radii of some of these neutron stars. For example Bogdanov & Grindlay (2009) found a lower limit R > 10.4 km (at 99.9% confidence) when Notes. The cameras and filters used, along with the total exposure time and the good time interval (GTI) used for the analysis are given in ks. The radii (source and background) for extracting spectra and light curves and the total counts in the source and the background regions are also provided.
considering a neutron star of 1.4 M for PSR J0030+0451. Published constraints on the radius for PSR J0437-4715 are R > 11.1 km (99.7% confidence) for a neutron star mass of 1.76 M (Bogdanov 2013), but more recent results suggest a mass of 1.44 ± 0.07 M (Reardon et al. 2016), which relaxes the lower limit slightly (Guillot et al. 2016). However, it is often difficult to determine accurate radii partly due to poorly constrained masses and the degeneracy between the mass and the radius. Precise knowledge of the mass improves the estimate of the radius, and well-constrained distance measurements also help somewhat (e.g. Lo et al. 2013). This in turn yields a more precise constraint on the neutron star equation of state. PSR J1909-3744 has been previously detected in the X-ray domain by Kargaltsev et al. (2012). The 29.7 ks Chandra ACIS observations revealed 63 net X-ray counts. The median energy of these photons is 1.0 keV, which could indicate that the emission from this pulsar is predominantly soft thermal emission. In order to determine whether this is indeed the case, we obtained XMM-Newton observations of PSR J1909-3744. It is a 2.95 ms pulsar which has both a well-constrained mass of 1.54 ± 0.03 M (1σ error, Desvignes et al. 2016) and a very well-constrained distance of 1.07 ± 0.04 kpc (Jones et al. 2017) thanks to its proximity and edge-on orientation allowing measurement of the range and shape of the Shapiro delay (Jacoby et al. 2005).
In this paper we investigate the nature of the X-ray emission with an aim to constraining the neutron star equation of state. Identifying new, thermally emitting pulsating millisecond pulsars with excellent mass and distance constraints will provide new targets for the NICER mission (e.g. Arzoumanian et al. 2014;Özel et al. 2016b), for XMM-Newton (Jansen et al. 2001), as well as for future missions, such as Athena (Nandra et al. 2013) or Strobe-X (Ray et al. 2019) all towards finally understanding the nature of the material inside neutron stars.

Data reduction and analysis
PSR J1909-3744 was observed with XMM-Newton on 20 September 2016 for a total of 51.6 ks. All three EPIC cameras were operated in full-frame mode. None of these observations have a high enough time resolution to detect any X-ray pulsations from this pulsar. Further information on the observations are given in Table 1. The optical monitor was not in operation during these observations.
We reduced the raw XMM-Newton data using the XMM-Newton Science Analysis System (SAS, version 16.1) and the latest calibration files at the time of the data reduction (CCFs, August 2017). The MOS data were reduced using the SAS task "emproc" and the SAS task "barycen" was used to barycentre the data, using the coordinates of the pulsar. The event lists were filtered with the #XMMEA_EM flag, and zero to 12 of the pre-defined patterns (single, double, triple, and quadruple pixel events) were retained. We identified periods of high background in the same way as described in the XMM-Newton SAS threads 1 and the good time interval is given in parentheses in Table 1. We also filtered in energy, using the range 0.2-12.0 keV. The pn data were reduced using the "epproc" and zero to 4 of the pre-defined patterns (single and double events) were retained, as these have the best energy calibration. Again we used the task "barycen" to barycentre the data. The background was treated in the same way as for the MOS data. We used the #XMMEA_EP filtering and the same energy range as for the MOS.
The SAS provides a task ("especget"), which allows the user to find an extraction region that optimises the source signal with respect to the background. We extracted the data using "especget" and the regions used are given in Table 1. The background was chosen from a source-free region close to the source and on the same CCD. These regions can be seen in Fig. 1. To create the spectra, we re-binned the data into 5eV bins as recommended in the SAS threads. We used the SAS tasks "rmfgen" and "arfgen" to generate a "re-distribution matrix file" and an "ancillary response file", for each spectrum. The pn data were binned to contain at least ten counts per bin and the MOS data to contain five counts per bin. The spectra were fitted using Xspec version 12.5 (Arnaud 1996).
The light curves were extracted using the same regions as for the spectra and using the maximum temporal resolution of the camera in use, as well as using the same temporal range for the source and background light curves for each camera. The source light curves were corrected for the background using the task "epiclccorr".

Results
PSR J1909-3744 was significantly detected with all three cameras, with 140 ± 12, 33 ± 6, and 38 ± 6 background subtracted counts detected with the pn, the MOS 1, and the MOS 2 respectively, see also Table 1. To determine if the X-ray emission from PSR J1909-3744 is thermal or non-thermal, we fitted simple   Top plot: data fitted with fixed absorption (tbabs) and a hydrogen atmosphere model (nsatmos, the * model in Table 2, green solid line). Bottom plot: same data divided by the model and shown as residuals, in units σ, which is the error bar of each data point. models to the MOS and pn spectra simultaneously. We chose a black-body model as a proxy for thermal emission and a powerlaw model to represent non-thermal emission. We also included a model for the absorption due to the interstellar medium, namely "tbabs" in Xspec with the Wilms et al. (2000) abundances. As the data were binned to contain only five to ten counts per bin, we used the C-statistic to assess the accuracy of the fit. Table 2 gives the results of the spectral fitting. The absorption is not constrained when using an absorbed black-body model. However, no absorption is not physical. Different H I maps indicate that we could expect n H ∼ 6.5 × 10 20 cm −2 (Kalberla et al. 2005) to ∼7.5 × 10 20 cm −2 (Dickey & Lockman 1990). Alternatively, using the conversion of dispersion measure to n H (He et al. 2013) and the dispersion measure of 10.4 pc cm −3 determined by Jones et al. (2017), we expect n H ∼ 3.1 × 10 20 cm −2 . We therefore used this value to fit the spectra and the model fits are also given in Table 2. We explore the degeneracy of the absorption and the temperature of the black body in Fig. 3 (top panel). As can be seen in Table 2, the power law with fixed absorption gave a poor fit (C-statistic of 52.5 as opposed to 31.4 for a black-body model with a fixed absorption component, see Table 2). However, the C-statistic does not allow us to directly determine the accuracy of the fit. We therefore did this in "Xspec" by simulating spectra based on the model 100 000 times, using parameter values drawn from a Gaussian distribution centered on the best fit. We then fitted each fake data set and calculated the test statistic. This goodness-of-fit testing only allows us to reject a model with a certain level of confidence (i.e., it can not provide a probability that the model is correct). Doing this for the fixed absorption and the black body, the model is rejected only at the 10.63% level, whereas the fixed absorption and the power law is rejected at 54.98%. This is not very conclusive, so in order to use χ 2 statistics we also tried fitting the data Fig. 4. 1, 2, and 3 σ confidence contours from fitting the X-ray spectrum of PSR J1909-3744. This shows the degeneracy between the temperature and size of the emitting surface when using the black-body model and the nsatmos model. For the latter, the model with fixed neutron star radius (R = 10.6 km, model (z) in Table 2) was used. Lines of constant (un-absorbed) bolometric flux are shown in grey, ranging from 10 −15 to 10 −13 erg s −1 cm −2 , in steps of 0.5 dex.
that was binned to contain a minimum of 20 counts per bin. The derived parameters we obtained when fitting the binned data are very similar to those determined when using the fixed absorption and either the black-body or the power-law models. For the fixed absorption and the black-body model we find a χ 2 ν = 0.67, with seven degrees of freedom and a null hypothesis that the observed data are drawn from the model of 0.7. For the fixed absorption and the power-law model we find a χ 2 ν = 1.98, with seven degrees of freedom and a null hypothesis that the observed data are drawn from the model of 0.05. This further reinforces the notion that the data are better fitted using a black-body model and that the emission is therefore predominantly thermal, with no need for a harder power-law tail as is sometimes seen in pulsar spectra, especially with low signal-to-noise ratio. Further, as can be seen in Fig. 3 (bottom panel), fitting the data with a power law requires a high photon index (reminiscent of a thermal spectrum) at the same time as a high absorption, much higher than the expected absorption. Both indicate that the power-law model is unlikely to describe the data.
We estimate the radius of the emitting region in a crude way by taking the emission over the whole observation, in order to have a spectrum with enough of a signal-to-noise ratio to fit with a model, and therefore smear out the effects of the projected spot area with rotation phase. Without the detection of pulsations, we also do not know if the emission is coming from one or two spots, so the emission radius is only an estimate and does not take into account that the spots are likely to be fairly flat, in addition to other factors outlined in Szary et al. (2017). We use the bbodyrad model in Xspec to estimate an approximate emitting radius of ∼58 m (see Table 2) typical of emission from polar caps with low signal-to-noise spectra (e.g. Bogdanov et al. 2006, and see Sect. 4). The emission is therefore coming from at least one polar cap. Given the low signal-to-noise, the degeneracy between the emitting area and the temperature of the emission is given in Fig. 4 for the black-body model. The same plot for the neutron star atmosphere model nsatmos is also shown in that same figure. Contours in Figs. 3 and 4 were obtained by sampling the parameter spaces via Markov Chain Monte Carlo in Xspec (command chain with the Goodman-Weare algorithm, 200 walkers, 100 000 steps, and a 20% burn-in). The un-absorbed luminosity of this MSP is 1.5 × 10 30 erg s −1 (0.2-10.0 keV), similar to other thermally emitting MSPs.
Fitting the X-ray spectrum with a neutron star atmosphere model (nsatmos in Xspec, Heinke et al. 2006), we obtain an equally accurate fit, as with the black body when we use the known distance and mass of the neutron star given in Sect. 1 (see Table 2). The data fitted with this spectral model can be seen in Fig. 2.
Whilst the best time resolution (73.4 ms for the pn in fullframe mode) is too coarse to detect any X-ray pulsations of the pulsar (2.95 ms), and the length of the observation, ∼9 h, is also too short to detect the binary orbital period of 1.5 d (Jacoby et al. 2005), we checked the background subtracted light curve for variability using a Kolmogorov-Smirnov test and a χ 2 probability of constancy test. Both tests gave a probability of constancy of between 0.1 and 0.01 (depending on the camera and the test), confirming that the X-ray light curves show no evidence for variability.

Discussion
Analysis of the X-ray spectrum shows that the X-ray emission from PSR J1909-3744 appears to be predominantly thermal. This millisecond pulsar has a very tight mass constraint, similar to that of PSR J1614-2230 which has a mass contraint of 1.928 ± 0.017 M (1σ, Fonseca et al. 2016). It also exhibits predominantly soft X-ray emission, making it a useful candidate for future studies to constrain the neutron star equation of state, especially if X-ray pulsations can be detected (see Sect. 1). However, the current short observations result in signal-to-noise ratios too low to make any such constraints.
The radius estimate of the X-ray emitting area is small compared to the size of a neutron star, which typically has a radius of 10-14 km. The emission is therefore likely to be coming from the polar caps of PSR J1909-3744, as is expected from such old objects where the surface of the neutron star has cooled so that it is no longer detectable in the X-ray band. The emitting radius is somewhat smaller than the classical radius of a polar cap, R pc = 2πR cP R (e.g. Dermer & Sturner 1994), where R is the neutron star radius, c the speed of light in a vacuum and P the rotation period of the neutron star. The latter is ∼2.9 km if we suppose a 10.6 km radius for the PSR J1909-3744 neutron star, the average of the best radius values for a neutron star of M = 1.5 M as determined by Özel et al. (2016b) who undertook a comprehensive study of 12 neutron stars. Considering the range of typical radii of neutron stars (10-14 km), the classical polar cap radius is between 2.7 and 4.4 km. Bogdanov et al. (2006) state that the radius determined through fitting the pulsar spectrum can be smaller than expected when the spectrum has been fitted with a single-temperature model, but that a higher signalto-noise spectrum would require a two-temperature model, as for PSR J0030+0451, PSR J2124-3358 (e.g. Bogdanov et al. 2008), and PSR J0437-4715 (Bogdanov 2013). Fitting the spectrum with two black bodies does not improve the fit, again probably due to the low signal-to-noise ratio. Further, it is known that the radii of neutron stars are underestimated when using a blackbody model as opposed to a realistic neutron star atmosphere model (e.g. Heinke et al. 2006). However, assuming the classical radius of the polar cap for this neutron star (2.9 km) in the nsatmos model provides a poor fit to the data with a C-statistic = 171.6, 33 degrees of freedom, suggesting that such a radius is too large. In fact, the nsatmos model produces a polar cap with size ∼100-400 m, slightly larger than the ∼50-300 m radius obtained with the black-body model (see Fig. 4).
The energy loss rate due to spin down corrected for the Shklovskii effect,Ė, is 4.3 × 10 33 erg s −1 (Smith et al. 2017). This value is similar to the values ofĖ obtained for the majority of the thermally emitting MSPs in 47 Tuc (Bogdanov et al. 2006) and in other thermally emitting MSPs (Abdo et al. 2013), which have also been detected in gamma-rays. It is therefore not surprising that this pulsar has just joined the list of Fermi LAT detected pulsars (Smith et al. 2017). Two gamma-ray peaks are observed, along with a very narrow radio peak (Smith et al. 2017). The X-ray luminosity (L X ) toĖ ratio for PSR J1909-3744 is similar to other thermally emitting millisecond pulsars, 4.4 × 10 −4 . The 47 Tuc pulsars have ratios between ∼10 −4 and 10 −3 (Bogdanov et al. 2006) and the non-globular cluster, thermally emitting pulsars have values of ∼10 −4 Abdo et al. (2013), Marelli et al. (2011). Modelling of the gamma-ray and radio light curves should eventually reveal further constraints on the pulsar geometry (e.g. Johnson et al. 2014).
The magnetic field at the light cylinder (r c ) can also be calculated if we assume a simple dipole model for the neutron star, where r c = cP/2π. We determine 7.4 × 10 4 G, which is again similar to values determined for the thermally emitting MSPs in 47 Tuc (Bogdanov et al. 2006) and average for Fermi LAT detected pulsars (Abdo et al. 2013).

Concluding remarks
To obtain an accurate constraint on the radius of the neutron star PSR J1909-3744, longer and higher time resolution observations are required to search for X-ray pulsations. If such pulsations are detected, the excellent mass and distance constraints combined with modelling of the X-ray light curve and spectra will allow strong constraints to be made on the neutron star equation of state. Whilst this source is likely to be too faint for NICER given that it is a factor of two fainter than PSR J1614-2230 and that there are sources at least as bright as PSR J1909-3744 at only 1 which would cause a high background for NICER observations, very long observations with XMM-Newton or shorter Athena observations will be able to achieve good enough quality spectra and light curves to constrain the radius to a few percent.