Redshifted absorption lines from an isolated neutron star

The solid crust constituting the outer layers of a hot neutron star is wrapped by an mm-to-cm thin atmosphere. Even if the atmosphere is so thin, it substantially affects the blackbody spectrum emitted by the surface, resulting in an overall hardening of the emitted spectrum 1 . The composition of the atmosphere has so far remained elusive. Several narrow absorption features have been detected and interpreted as arising from proton (or electron) resonant cyclotron absorption in the neutron star magnetic ﬁeld 2–4 . Apart from these, for a Hydrogen atmosphere no spectral features are expected, whereas when it is polluted with metals, absorption features start appearing in soft X–ray spectra 5,6 . Absorption edges and features have been possibly observed during thermonuclear explo-sions onto the neutron star surface 7,8 . Isolated neutron stars represent a breeding ground where to look for absorption features, thanks to their simple X–ray spectra. Here we report on the detection of redshifted Nitrogen and Oxygen absorption features from the closest and brightest isolated neutron star. The lines are ∼ 50 eV wide and their intensity is incompatible from originating in the interstellar path to the neutron star. Lines are redshifted by a common gravitational redshift of z g = 0 . 216 ± 0 . 004 . The magniﬁcent seven is a group of nearby, cooling, relatively young, and isolated neutron stars 3,9 . Their X–ray spectra are very simple, consisting of a cooling thermal component.

struments, namely the Reflection Grating Spectrometers (RGS 17 ). After a standard filter and data extraction procedure (see Methods for more details), we obtained 2.31 Ms and 2.24 Ms exposures for the RGS1 and RGS2, respectively, comprising 0.38 Mcounts in the 0.324-0.8 keV energy range for the RGS1, and 0.33 Mcounts in the 0.332-0.5 plus 0.625-0.8 keV energy ranges (due to the failure of one CCD) for the RGS2. To complement these data, we also added Chandra Low Energy Transmission Grating Spectrometer (LETG 18 ). We collected 0.61 Ms with the High Resolution Camera (HRC,0.18 Mcounts in the 0.22-0.28 keV plus 0.32-0.8 keV energy range) and 0.25 Ms with the Advanced CCD Imaging Spectrometer (ACIS, 6 kcounts in the 0.22-0.28 keV plus 0.41-0.8 keV energy range). Given our interest in narrow absorption (or emission in principle) features, we do not focus too much on the overall continuum model: we just need to have a good fit of the continuum on top of which searching for lines. This has been achieved with an absorbed double blackbody model (a blackbody and a neutron star atmosphere model provides similar results). Given the superb statistics of our spectra, we leave the column density, blackbody temperatures and radii for each instrument to vary independently. We obtained broadly consistent results for each instrument (see Methods).
A first iteration is then applied on the ungrouped spectra to get rid of small-width (< 1 eV, Gaussian sigma) absorption/emission lines, using Cash statistics to remove possible spurious features. We found a few of them. We then turned to χ 2 statistics grouping data to 10 (RGS1, RGS2), 300 (LETG-HRC), and 150 (LETG-ACIS) counts per energy bin (the higher grouping for the Chandra data is driven by the higher background and spectral resolution) and Churazov weighting scheme to use χ 2 statistics. We then added iteratively Gaussian lines (either with a positive or negative normalisation in the 0.25-0.8 keV energy range) until the variation in the χ 2 is more than 10. This turned out to be a very good prescription to identify high significance lines (see Methods for more details).
We end up with two lines. The first line is at 411 ± 1 eV (errors are computed for ∆χ 2 = 2.7, equivalent to a 90% confidence level for one parameter of interest) and its inclusion in the fit produces a ∆χ 2 = 46.2. The line is structured and it is modelled with two Gaussians, one narrow 4±2 eV wide and the other large 38±10 eV wide. The equivalent width of the line is 4.2±2.3 eV (3.5 eV and 0.6 eV for the large and narrow components, respectively). The second line is at 562±10 eV (∆χ 2 = 30.2), it is 32±12 eV wide, and its equivalent width is 6.5±2.2 eV.
The identification of the two lines is not straightforward. The low energy line can be tentatively identified with NV Lyα (also known as Kα, at 421.5 eV), only including a non-negligible redshift of ∼ 0.03, and the high energy line with OVI Lyα (at 563.1 eV) with no redshift (see Fig. 1). Line energy is not the only important quantity, however. Converting the line equivalent width to a, more physical, column density with a curve of growth analysis 19, 20 , one obtains very large values (N OVI ∼ 5 × 10 17 cm −2 and N NV ∼ 3 × 10 17 cm −2 ) for each element. This value has been computed entirely attributing the line width to the thermal broadening (∼ 20, 000 km s −1 ), thus providing the lowest possible column density. Even assuming that these are the only atomic species, assuming Solar abundances 21 , one obtains equivalent column densities of N H ∼ 10 21 cm −2 and N H ∼ 4 × 10 21 cm −2 , respectively. These values are much larger than the observed equivalent column density along the line of sight to RXJ1856, which is a few 10 20 cm −2 . This is true for any kind of material, diluted or compact, along the line of sight.
Weak absorption features have been observed in all the seven isolated neutron star, except of RX J1856: narrow features around 0.57 keV and 0.53 keV were observed in RX J0720.4-3125 4, 23 the second brightest isolated neutron star; RX J1605.3+3249 showed a broad absorp-tion feature at 0.40 keV and a narrow feature at 0.57 keV 24-26 ; RX J1308.6+2127 showed an absorption feature at 0.53 keV in the co-added RGS spectra 4 . Phase dependent absorption features were observed in RX J0720.4-3125 and RX J1308.6+2127, and interpreted as proton cyclotron lines 27 . Given the magnetic field of RXJ1856 (∼ 2 × 10 13 G), this feature is expected to lie outside the relevant energy range (∼ 0.06 keV for protons and ∼ 120 keV for electrons).
Using a model of a condensed iron surface, with a partially ionised thin hydrogen atmosphere a gravitational redshift of z = 0.15 and z = 0.205 has been indirectly estimated in RBS 1223 and RX J0720.4-3125, respectively 28, 29 .
The detection of absorption features close in energy to the ones reported here testifies that no calibration issues are present. These lines were usually interpreted either as arising from the hot interstellar medium or in a localised hot medium, possibly clumpy, along the line of sight. The observations of similar features in a number of isolated neutron stars make the second hypothesis untenable: the hot medium should have a large covering factor increasing too much the overall mass 4, 23 . The close distance to and the lack of detected features in RXJ1856 was indeed taken as a proof for the hot interstellar medium hypothesis 4 . Our line detection in the closest isolated neutron star (∼ 120 pc) indicate instead that also this scenario has to be discarded, leaving only lines originating at the star surface and redshifted by gravity.
Following this line of reasoning, we are forced to attribute these lines to the atmosphere of the neutron star, where matter is not lacking, and the emission spectrum is produced. In this case, lines will be redshifted by the neutron star gravitational field and also the atmosphere temperature should be higher, allowing only for the highest ionisation species to be present.
Thus, our line identification has to be changed. We consider the two H-like ions NVII Lyα (at 500.3 eV) and OVIII Lyα (at 653.6 eV). We fix the line energy and allow for the same (free) redshift of the lines. We also require that the two wide lines have the same width. We derived a comparable fit and a gravitational redshift z g = 0.216 ± 0.004. The Gaussian line with is 54.6 +6.8 −5.7 eV. Considering the full width at a half maximum we obtain ∆λ/λ ∼ 0.25.
At this point we searched for a confirmation in the XMM-Newton pn data. As already stated, pn data suffer from calibration issues (when challenged with this huge number of counts, almost 6 millions) and the spectra are position dependent. Following Sartore et al. 16 , we identified five positions on the detector, extracted the relevant spectra and response files (see Methods). We then fit the data with an absorbed double blackbody model, finding four weak and broad absorption features. Two of them are consistent with those found in high resolution spectral data. The equivalent widths are still too large (∼ 10 eV) to come from material along the line of sight. Looking for redshifted lines we tentatively identify NVII Lyα, NVII Lyβ, OVIII Lyα, and OVIII Lyβ, at a common redshift of z g = 0.224 ± 0.015. This value is consistent with the one derived from high resolution spectral data. Given the uncertainties in the positiondependent calibration of the pn detector, however, we consider this as an indirect confirmation of the results obtained with high resolution spectral data.
The value of z g calls for standard values for the neutron star paramaters: for a standard mass of 1.4 M we derive a radius of 12.8 ± 0.2 km or, for a 12 km radius a mass of 1.31 ± 0.01 M (Fig. 2). This might indicate that magnetars are born with a mass somewhat lower than the canonical neutron star mass. Atmosphere models with a mixture of elements, as foreseen by our spectral fits, do not exist and should be developed to test in more depth our results. This can lead to derive directly from the fit the mass and radius of the neutron star. A 100 ks Athena X-IFU observation will provide a ∼ 1% error on the gravitational redshift and a comparable precision on the independent estimate of the neutron star mass and radius. The constraint from the two redshifted absorption lines is shown with a light blue solid region.

Methods
Data preparation We collected 42 XMM-Newton observations (see Table S1). We processed all of them following the analysis threads on the XMM-Newton Standard Analysis Software (SAS) pages 1 . We used SAS 18.0.0 and the latest available calibrations. We filtered the soft proton flare background using a threshold on (standard) rate of 0.2 counts s −1 , computed on the RGS1 data. Background and response matrices were also generated. All spectra were then visually inspected and summed on an instrument basis, together with background and response, using rgscombine. We end up with 2.31 Ms (2.24 Ms) for the RGS1 (RGS2), summing 383,021 (325,078) counts in the 0.324-0.8 keV (0.332-0.5 plus 0.625-0.8 keV, due to the failure of one CCD) energy range. Second order spectra are background-dominated (∼ 8% of source photons) and were discarded. We decided to limit the spectra to 0.8 keV not to encounter problems with the background and because beyond this energy, the number of photons starts becoming too low to search for narrow features. Data were binned to 10 counts per energy bin, given the large number of photons and the relatively low background (82%). A Churazov weighting scheme has been adopted in XSPEC (12.10.1f) to ensure χ 2 statistics 32 .
Simultaneously, we processed the pn data taking the (standard) threshold on rate of 0.4 counts s −1 , to filter out soft background flares (being the pn more sensitive than the MOS CCDs working in the RGS, we adopted this more restrictive filtering). All pn observations were collected in small window mode, apart from 4 observations during which the pn was not working or in timing mode. Spectra were extracted from a 880 pixel radius centred on source. Background spectra were extracted from a 600 pixel circular region at the edge of the small window. Single and double events were retained and FLAG==0 was applied. Ancillary  and response files were generated using arfgen and rmfgen and the latest calibration files.
Following Sartore et al. 16 , we sum pn spectra on a detector position basis. We individuated five groups of positions (see Table S1). We collected 5.1 Mcounts in 2.0 Ms in the 0.2-0.8 keV energy range. Data were binned to 25 counts per energy bin. A few observations did not find a suitable group to join (thus having a small statistics) and were discarded.
For the Chandra data, we used CIAO v.4.12 and CALDB 4.9.1. We reprocess all the public data with chandra repro and follow the analysis threads 2 . We collected data for 0.61 Ms for the HRC (6 observations) and 0.25 Ms for the ACIS (10 observations, see Table   S2). Only first order data (plus and minus) were considered and summed up. We restricted the energy range of the LETG-HRC to 0.22-0.28 keV plus 0.32-0.8 keV to get rid of the effective area decrease around the carbon edge. The total number of collected counts is 0.18 Mcounts.
For the LETG-HRC the background (68% of the total) and the spectral resolution are high, so that we binned the data to 300 counts per energy bin. LETG-ACIS data can be used only when the ACIS-S is on the grating path, due to its higher sensitivity to soft photons. We restricted the energy range to 0.22-0.28 keV and 0.41-0.8 keV, where we collected 6 kcounts, mainly in the very soft part of the spectrum. Data were binned to 150 counts per energy bin.
Before starting with the line fitting procedure, we fit the single spectra binned to one photon per energy channel (and adopting the Cash statistics) with a double blackbody model.
We then add a Gaussian single line with width (σ) fixed to 0.5 eV. This line was then left free to vary across the entire spectrum (0.2-0.8 keV) to get rid of single channels, coincident with hot pixels or mis-calibrated channels.
Line search: observer-frame lines After grouping and filtering we started with the search for a meaningful continuum model over which search for spectral features. Spectral fits were carried out with the spectral package XSPEC 30 . After some attempts, we considered a model made by two blackbody spectra (one soft and one harder) absorbed by a column density. The fit is equivalent to the one that we can obtain from an absorbed blackbody plus a neutron star atmosphere model, but it is faster in the computation (being analytical and not relying on interpolations). The absorbing column density has been modelled with tbabs 21 , using Verner cross-sections 31 and wilm abundances. Given the large number of counts and small intercalibration uncertainties, we leave free all the parameters for each instrument. The blackbody radius has been limited to 20 km and the distance fixed to 120 pc. Only for the fit of LETG-ACIS data the soft blackbody component is unconstrained. The results of the fit are reported in  We also tried the more complex tbnew absorption model. We tried to leave free all the element abundances, but only N and O were significantly different from 1. Leaving free the abundances of these two elements, we obtained an improvement in the fit of ∆χ 2 = 17.5 (and a χ 2 r = 1.22 for 946 degrees of freedom . These abundances are much higher than solar and are not easily interpretable. In addition, the fit is still not statistically acceptable. We then started the iteration process adding a Gaussian line (using tbabs to model interstellar absorption). The line was kept free to vary in the 0.25-0.8 keV energy range, with a free Gaussian width in the 1-100 eV range. We decided to stop when the the improvement in the χ 2 for the addition of a new line is less than 10. After a first fit, a search of the best fit line energy is obtained performing a steppar search with 1199 or 1201 steps (i.e. ∼ 0.5 eV) resolution. The variation of the χ 2 is annotated (see Table S4) and all the line parameters are left free to vary when we include a further line. As a visual impression, in Fig. S1 we show the improvement in the χ 2 as a function of energy depending on three fixed line widths (in the fit the line width is left free to vary).
We also report in Table S4 the χ 2 variation on an instrument basis. We report in Table   S5 the final line values. The addition of a further line produces a ∆χ 2 = 2.3 and we deemed it a statistical fluctuation. This putative line has an energy close to the centroid of Line2 and has a line width that saturates our limit of 100 eV. It is apparent from the values in the upper part of Table S4 that Line1 and Line3 are very close (and consistent in energy) and that RGS1 and RGS2 react in different ways to the inclusion of this line. We therefore restart the process different from zero at > 3.5 σ confidence level. The overall fit is reported in  Line search: redshifted lines We searched for redshifted lines. We considered NVII Lyα at 500.3 eV and OVIII Lyα at 653.6 eV as the two candidates. These lines were identified following an iterative process with the main lines (C, N, O and Ne) to identify which can be coupled to a common redshift. We then fit the data using a redshifted Gaussian zgaus instead of a Gaussian, keeping the structured line for NVII Lyα. In addition, we impose the same line width to the two wide lines. Given the smaller degrees of freedom we expect a slightly worse fit with respect to the unredshifted lines. We obtain a χ 2 = 1103.8 for 943 dof (χ 2 r = 1.17).
Leaving free the two line widths improves the χ 2 by 6.4, resulting in a probability of chance coincidence by means of an F-test of 11.6%. If we instead leave the two redshifts independent, the χ 2 improves by 9.7, with an F-test probability of 1.9%. Fit parameters are reported in Table S6.   Table S5: Fit with an absorbed double blackbody model and lines. Errors were computed with ∆χ 2 = 2.71, or 90% confidence level for one parameter of interest. The normalisation K is related to the equivalent emitting radius R by the relation K = (R/10km) 2 /(d/10kpc) 2 , with d the source distance, 120 pc. A * close to one parameter means that it is consistent with zero. The equivalent widths refer to RGS1. and N NVII ∼ 2 × 10 20 cm −2 . If distributed isotropically in the neutron star atmosphere this give ∼ 2 × 10 11 g of OVIII and ∼ 10 11 g of NVII, assuming 12 km radius and 1 cm thickness.
While there is a desert of strong high-ionisation lines around the energy of the N VII Lyα, this is not the case for the OVIII Lyα. One can in principle envisage that the line identification is not correct and this second line might be attributed to OVII Lyβ (at 665.6 eV), or even to OVII Lyγ (at 697.8 eV). We first fix the line energy of the second line to OVII Lyβ (at 665.6 eV). The fit improves to χ 2 = 1095.0 for a ∆χ 2 = 8.9. Based on the Akaike information criterion (AIC = 2k + χ 2 , where k is the number of parameters 33 ), the OVII Lyβ improves the fit at a significance level of (equivalent) 2.5 Gaussian σ. We also left free the energy of this second line as well as the redshift and anchor the fit to the N VII Lyα. In this case the line energy ends at 667 +16 −11 eV, in between OVII Lyβ and OVII Lyγ. Based on the AIC criterion, the fit improves by 2.6 Gaussian σ. Independently on the line identification, as a result of these experiments the gravitational redshift move to 0.217 ± 0.004, with a mere 0.5% change.
Line search: pn detector As shown by Sartore et al. 16 , the spectrum of RXJ1856 on the pn detector is position dependent. This is mainly due to the huge number of counts coupled to small calibration uncertainties. We identified five groups of positions (see Table S1, A to E).
Data were grouped on a position basis and spectra were independently fitted with the absorbed double blackbody model, adding 0.5% systematic error and allowing for a 20 km neutron star maximum radius. For four of the spectra the second blackbody component was not needed (group C needed it). The column densities range from (0.3 ± 0.1) to (0.9 ± 0.6) × 10 20 cm −2 and the blackbody temperature 61.3 ± 0.1 to 63.4 ± 0.1 eV. The overall spectrum provides a reduced χ 2 r = 1.68 for 552 dof. Residuals are shown in Fig. S3. We iteratively add lines, as   The variation in the χ 2 is ∆χ 2 = 54.5. The forth line is at 570.2 ± 7.1 eV and the equivalent width of 17.4 ± 1.8 eV. The variation in the χ 2 is ∆χ 2 = 25.6. The overall fit is presented in Fig. S4 and line energies are reported in Table S7. The lines at 411.5 eV and 570.2 eV can be identified with the lines found in the high resolution spectral data. The other two lines found in the on analysis did not pass our threshold during the analysis of high resolution data. We go back to the RGS and LETG data to derive an upper limit for the lines at 661.1 eV and 490.1 eV. The inclusion of the first (fixed energy) line provides a ∆χ 2 = 7.1 and the upper limit on the equivalent width is < 12.9 eV. The inclusion of the second line provides a ∆χ 2 = 4.6 and the upper limit on the equivalent width is < 14.8 eV. Both equivalent widths are slightly smaller than the values derived from the pn data.
As before, these lines can be tentatively identified as NV Lyα (412 ± 4 eV vs. 421.5 eV), NVII Lyα (490 ± 5 eV vs. 500.3 eV), OVII Lyα (570 ± 7 eV vs. 574.0 eV), and OVIII Lyβ (661 ± 9 eV vs. 665.6 eV), respectively, even if the dependence on position of the pn spectral response, make us not fully trust the pn spectra. As discussed before, however, the line equivalent widths are too intense to come from the interstellar medium and the Nitrogen lines require the same (small, z ∼ 0.02) redshift, which is not requested by Oxygen lines. As above, we are forced to look for gravitationally redshifted lines.
Given the high temperature of the neutron star atmosphere, we were driven in our choice by the lowest level of the highest ionisation level species. Tentatively, we identify the lines with NVII Lyα (the line with observed energy at 412 eV), OVII Lyα (490 eV), OVIII Lyα (570 eV), and OVIII Lyβ (662 eV). The overall fit is good with a reduced χ 2 = 0.97 with 544 dof (67% nhp). The line width is 46 ± 5 eV. The equivalent widths are 18.7 ± 10.6 eV (NVII Lyα), 9.0 ± 4.2 eV (NVII Lyβ), 16.8 ± 4.1 eV (OVIII Lyα), and 13.7 ± 8.2 eV (OVIII Lyβ). The gravitational redshift is z g = 0.224 ± 0.015. This value is consistent with the one derived from high resolution spectral data. Given the uncertainties in the position-dependent calibration of the pn detector, we consider this as an indirect confirmation of the results obtained with high resolution spectral data.  Given that data were collected over several years, the timing solution is not good enough to afford a pulse phase spectroscopy study.
MOS data were all taken in timing mode and the usable energy range starts from ∼ 0.5 keV. For this reason MOS data were not considered here.
Line search: tentative physical modelling We also tried a physical modelling of the absorber using XSTAR. XSTAR is used for calculating the physical conditions and emission/absorption spectra of photoionised gases 34 . We used the analytical version of XSTAR, i.e. warmabs (v.2.2, 35 ). We adopted the highest population for a density of 10 16 g cm −3 . Modelling the neutron star atmosphere is beyond the scope of this paper, so we approached the fit with this model having clear in mind that this cannot provide a completely satisfactory description of the absorber. After a first fitting run, it became clear that the hot column density is tightly coupled with the Oxygen and Nitrogen abundances, so we fixed it to its highest value, 10 24 cm −2 .
The overall fit is better than the one with an absorbed double blackbody model but worse than the one obtained with the addition of three Gaussian lines, as expected being the grid not able to perfectly model the large line widths, even if leaving free the turbulent velocity parameter improves the fit considerably. The fit values are reported in Table S8. The χ 2 is 1128.1 for 943 dof (χ 2 r = 1.20). The logarithm of the ionisation parameter ξ = L/(n r 2 ) (where L is the ionising luminosity, r is the distance, and n is the particle density) is high with a lower limit of 1.6 (90% confidence). The abundances of all elements but Nitrogen and Redshift (z) 0.214 ± 0.004 Table S8: Fit with an absorbed double blackbody model and hot absorbing medium. Errors were computed with ∆χ 2 = 2.71, or 90% confidence level for one parameter of interest.
A * indicates a parameter whose error is not constrained: we reported the value of the fit and the lower limit in parenthesis.