© ESO 2018

1 Introduction

X-rays have been detected from more than one hundred rotation-powered pulsars. The properties of the X-ray emission, and in particular the relative strength of the non-thermal and thermal spectral components, depend on the pulsar’s age (Becker 2009; Harding 2013).

Young pulsars have the highest surface temperatures, but their thermal flux is unobservable; they are outshined by a much stronger non-thermal emission, which fades as the pulsars age and slow down, thus unmasking the underlying thermal emission of the neutron star (NS) surface.

Middle-aged pulsars show a mix of non-thermal and thermal X-ray emission. The latter is often modelled with two blackbody components. The cooler components have emitting areas consistent with emission from the whole NS surface and temperatures in agreement with the prediction of standard cooling theories. The hotter components have blackbody temperatures typically in the range 0.1–0.3 keV and emitting radii between a few tens and a few hundred metres, implying that the X-rays come from a small region of the NS (i.e. PSR B0656+14, PSR J0633+1746, PSR B1055−52 (De Luca et al. 2005) and PSR J1740+1000; Kargaltsev et al. 2012).

X-rays have been detected also in older pulsars, that have cooled down to temperatures too low to significantly emit in the X-ray band. In these objects both non-thermal X-rays of magnetospheric origin and/or emission from hot polar caps have been observed with high confidence only in a few objects (see for example PSR B1929+10; Misanovic et al. 2008). In fact, most old pulsars are rather faint X-ray sources and their study can benefit of the use of maximum likelihood (ML) techniques, as demonstrated for PSR B0943+10 (Hermsen et al. 2013; Mereghetti et al. 2016) and PSR B1822−09 (Hermsen et al. 2017).

Here, we report a reanalysis of XMM–Newton archival data of four pulsars with characteristic ages in the range 0.3− 5 Myr based on ML methods. The main properties of our targets are summarized in Table 1.

2 Observations and data reduction

We used the data obtained with the pn (Strüder et al. 2001) and the two MOS (Turner et al. 2001) cameras, which constitute the European Photon Imaging Cameras (EPIC) instrument on board XMM–Newton and cover the energy range 0.2–12 keV.

The medium optical filter was used for PSR B0114+58 and for the MOS cameras data of PSR B0628−28. All the other data were obtained using the thin filter. The time resolution of the different operating modes used during our observations are indicated in Table 2. We note that while the pn time resolution was adequate to reveal the pulsations in all our targets, this is not true for the MOS data, except in the case of PSR B1133+16. The data reduction and analysis were done using version 15 of the Science Analysis System (SAS)¹ and XSPEC(ver. 12.8.2) for the spectral fits.

We reprocessed the pn data using the SAS task epreject to reduce the detector noise at the lowest energies. We then removed the time intervals of high background by rejecting all the periods with a pn count rate higher than 1.0 cts s⁻¹ in the range 10–12 keV. The resulting net exposure times for each pulsar are indicated in Table 2.

We used single- and multiple-pixel events for both the pn and MOS. The events detected in the two MOS were combined into a single data set, and in their analysis we used averaged exposure maps and response files.

We applied a maximum likelihood (ML) analysis to extract the net (i.e. background-subtracted) source counts used for the spectral andtiming analysis. The resulting count spectra were put in a format compatible with analysis in XSPEC. By exploiting the knowledge of the instrumental point spread function (PSF), this method provides spectra and light curves with a number of source counts larger than what is typically obtained with the traditional method based on “source” and “background” extraction regions. Another advantage is that the background is directly measured at the position of the source, contrary to what is done with the traditional analysis in which it is estimated from “source free” regions of the image.

The ML method measures the source counts taking into account all the events of the region of interest. Specifically, it finds the source flux and background level which maximize the probability of obtaining the number of total counts N_ij measured in each spatial pixel. The values of N_ij follow a Poisson probability distribution, with expectation value μ_ij = β + σ ×PSF_ij, where β is the background level, which we assumed to be uniform; σ is the number of source counts; and PSF_ij gives the relative contribution of the source flux in the pixel (i, j) taking into account the instrumental PSF, which depends on the position in the field of view and on the energy of the event. We used spatial pixels of 1 × 1 arcsec² and we analysed separately the data of the pn and of the sum of the two MOS, adopting the respective PSFs, derived from in-flight calibrations², with parameters appropriate for the average energy value in the considered bin.

Following Hermsen et al. (2017), we generalized the above method to also take into account the pulse phase information of the events. In this 3D-ML approach, the events are binned in spatial and phase coordinates and the expectation value of bin (i, j, k) is μ_ijk = β + σ_u ×PSF_ij + σ_p ×PSF_ij ×Φ_k.

Now σ_u and σ_p represent the source counts for the unpulsed and pulsed components, while Φ_k is the normalized pulse profile at phase bin k. With this approach, and assuming a sinusoidal pulse (consistent with the data of our targets), we could simultaneously derive the fluxes (and spectra) for the unpulsed and the pulsed components. We define the pulsed fraction as the ratio of the pulsed to the total counts, σ_p ∕(σ_p + σ_u).

3 Results

3.1 PSR B0628−28 (J0630−2834)

We applied the ML analysis to a circular region of radius 60″ centred at RA = 06^h30^m51^s.01, Dec = −28°35′06. ″8. In the energy range 0.3–10 keV, this yielded 1021 ± 40 source counts in the pn and 527 ± 28 in the MOS, corresponding to a total count rate of (4.0 ± 0.1) × 10⁻² cts s⁻¹.

The spectral analysis was carried out in 20 energy bins for the pn and in 13 for the MOS. The two spectra were fitted simultaneously. Among single-component models, an absorbed power law gave the best fit with $\Gamma ={2.85}_{-0.15}^{+0.16}$ and N_H $=({8.5}_{-2.1}^{+2.3})\times {10}^{20}$ cm⁻² (${\chi }_{\nu }^2=0.67$ for 30 degrees of freedom (-d.o.f.)). The error regions of the parameters are shown in Fig. 1a. A fit with a single blackbody was not acceptable (${\chi }_{\nu }^2=3.9$ for 30 d.o.f.). The best fit power law spectrum, together with its residuals and those of the blackbody fit are shown in Fig. 2.

The best fit N_H is within 1σ of the value of 10²¹ cm⁻², derived from the pulsar dispersion measure with the usual assumption of a 10% ionization of the interstellar medium (He et al. 2013). Fixing N_H to this value, we get tighter constraints on the photon index, Γ = 2.95 ± 0.06. The corresponding unabsorbed flux in the 0.2–10 keV energy range is f_0.2−10 = (1.08 ± 0.05) × 10⁻¹³ erg s⁻¹ cm⁻².

All the above results are in agreement with those of standard spectral analysis of the same data. Becker et al. (2005) in fact found a power law with photon index $\Gamma ={2.63}_{-0.15}^{+0.22}$, while Tepedelenlioǧlu & Ögelman (2005), simultaneously fitting the XMM–Newton and Chandra spectra, obtained $\Gamma ={3.20}_{-0.23}^{+0.26}$ (1σ errors in both cases). These authors also considered a power law plus blackbody model, obtaining acceptable fits with slightly harder power laws and blackbody temperatures kT ~ 0.3 keV.

Using a power law plus blackbody model we found that several combinations of blackbody temperature and normalization, which give only a limited flux contribution, are consistent with the data. The best fit parameters were Γ = 2.92 ± 0.07 and kT ~ 0.02 keV, while the emitting radius is poorly constrained. However, the addition of a blackbody is not statistically required. Therefore, we derived an upper limit on its intensity as follows. We fitted the spectrum with a power law plus a blackbody model (N_H fixed to 10²¹ cm⁻²) and computed the 3σ confidence ranges of the blackbody normalization and temperature, leaving the power law parameters free to vary. The results are shown in Fig. 3; the upper panel refers to the emitting radius and the lower panel to the bolometric luminosity. The latter has been computed as L_bol = πR²σT⁴ because we are considering emission from a hot spot on the stellar surface.

To investigate the spectral properties of the pulsed and unpulsed emission, we performed a 3D-ML analysis on the pn data in seven energy bins in the range 0.3 keV10 keV. The phases of the detected counts were computed using the period at the epoch of the XMM–Newton observation, P = 1.2444225976(2) s, obtained from the ephemeris reported in the online³ ATNF Pulsar Catalogue (Manchester et al. 2005). The pulsed fraction is higher in the 0.40–0.65 keV range, with an average value of 0.46 ± 0.05 and decreases at lower and higher energies (Fig. 4, left). The energy dependence of the pulsed fraction is also reflected in the different shape of the pulsed and unpulsed spectra (Fig. 4, right). The unpulsed spectrum is fitted well by a power law with the same slope as that of the total spectrum, while that of the pulsed emission requires a significantly steeper power law or a blackbody (Table 3).

Fig. 1 Confidence regions (1, 2, and 3σ) of the best fit parameters for the spectra of (a) PSR B0628−28, (b) PSR B0919+06, (c) PSR B0114+58, and (d) PSR B1133+16. The best fit values are indicated by a plus and the DM-based value of NH by a dashed line.

Fig. 2 EPIC-pn (blue diamonds) and -MOS (red squares) X-ray spectra of PSR B0628−28. The best fit power law is shown in the top panel, and the corresponding residuals in the middle panel. The bottom panel shows the residuals obtained by fitting the spectra with an absorbed blackbody.

Fig. 3 Upper limits (3σ confidence level) on the blackbody component in PSR B0628−28 as a function of blackbody temperature. The top panel gives the limits on the emitting radius, while the lower panel gives the corresponding limits on the bolometric luminosity computed as Lbol =πR2σT4.

Fig. 4 Left: pulsed fraction of PSR B0628−28 as a function of energy obtained with the 3D-ML analysis. Right: EPIC-pn total (black circles), unpulsed (blue diamonds), and pulsed (red squares) spectra of PSR B0628−28.

3.2 PSR B0919+06 (J0922+0638)

The faint X-ray counterpart of PSR B0919+06 was first reported by Prinz & Becker (2015). Based on a spectrum with 171 ± 18 counts, they rejected the blackbody model and found a good power law fit with $\Gamma ={2.3}_{-0.4}^{+0.8}$. Our ML analysis of the same data was carried out in circular region with radius of 45″ centred at RA = 09^h22^m14^s.20, Dec = +06°38′33. ″7. In the energy range 0.2–10 keV we obtained 223 ± 25 (pn) and 91 ± 15 (MOS) source counts, corresponding to a total count rate of (1.3 ± 0.1) × 10⁻² cts s⁻¹.

We simultaneously fitted the pn (eight energy bins) and the MOS (seven energy bins) spectra. The single-component model that best fits the data is an absorbed power law (see Fig. 1b and Table 4). The single blackbody model is statistically unacceptable (${\chi }_{\nu }^2=1.92$ for 12 d.o.f.). We show the comparison between the power law and the blackbody best fits in Fig. 5. Fixing the absorption to the DM-based value 8 × 10²⁰ cm⁻², we obtain Γ =3.1 ± 0.2 and an unabsorbed flux f_0.2−10 = (3.0 ± 0.4) × 10⁻¹⁴ erg s⁻¹ cm⁻².

Although a good fit is also obtained with a power law plus blackbody (Table 4), the addition of a blackbody is not statistically required. Therefore, as we did for PSR B0628−28, we derived the upper limits on its intensity as a function of T, which is shown in Fig. 6.

Fig. 5 EPIC-pn (blue diamonds) and -MOS (red squares) X-ray spectra of PSR B0919+06. The best fit power law is shown in the top panel, and the corresponding residuals in the middle panel. The bottom panel shows the residuals obtained by fitting the spectra with an absorbed blackbody.

3.3 PSR B0114+58 (J0117+5914)

This faint pulsar, observed with XMM–Newton for less than 6 ks, is the one with the smallest number of counts in our sample. Prinz & Becker (2015) obtained a poorly constrained spectrum, which could be fitted equally well by a power law with Γ =3.3 ± 0.5 or by a blackbody with temperature kT = 0.18 ± 0.03 keV and emitting radius R = 322 ± 161 m (for d = 1.77 kpc).

We applied the ML analysis in a circular region centred at RA = 01^h17^m38^s.24, Dec = +59°15′07. ″8 with a radius of 60″, from which we excluded a circle centred in RA = 01^h17^m34^s.87, Dec = +59°15′15. ″4 with a radius of 20″ to avoid a nearby source. This resulted in 85 ± 13 pn counts and 44 ± 9 MOS counts in the range 0.2–10 keV, corresponding to a total count rate of (2.2 ± 0.3) × 10⁻² cts s⁻¹.

Simultaneously fitting the pn and the MOS spectra it was impossible to discriminate between a power law (${\chi }_{\nu }^2=0.58$ for 7 d.o.f.) and a blackbody (${\chi }_{\nu }^2=0.62$ for 7 d.o.f.). However, restricting the analysis only to the pn spectrum, which is of better statistical quality, and fixing the absorption to the DM-based value of 1.5 × 10²¹ cm⁻², we found that the blackbody is preferred (${\chi }_{\nu }^2=0.32$ for 4 d.o.f., compared to ${\chi }_{\nu }^2=1.80$ for the power law). The best fit temperature is kT = 0.17 ± 0.02 keV, the radius of the emitting area is $R={405}_{-90}^{+110}$ m, and the 0.2–10 keV unabsorbed flux is f_0.2−10 = (4.7 ± 0.7) × 10⁻¹⁴ erg s⁻¹ cm⁻² (see Table 5 and Fig. 1c). Leaving the absorption free to vary, the power law is also acceptable, but with an unphysical photon index $\Gamma ={6.5}_{-2.0}^{+3.0}$.

3.4 PSR B1133+16 (J1136+1551)

PSR B1133+16 has both optical (Zharikov et al. 2008; Zharikov & Mignani 2013) and X-ray (Kargaltsev et al. 2006) counterparts. According to these works, its X-ray spectrum can be described by either a power law or a blackbody (or a combination), but the statistics are too poor to distinguish between these models.

PSR B1133+16 was observed with five separate pointings over a time of about one month. With a simultaneous fit of the separate spectra, Szary et al. (2017) found a power law with Γ ≈ 2.5 and reported that the addition of a blackbody with $kT={0.25}_{-0.03}^{+0.05}$ keV and a radius of ${14}_{-5}^{+7}$ m (for d = 350 pc) is also consistent with the data.

We verified that the analysis of the single observations gave fully consistent results. Therefore, we report only those obtained by combining the five pointings to extract a single pn and a single MOS spectrum. We applied the ML analysis in a circular region of radius 35″ and centred at RA = 11^h36^m02^s.527, Dec = 15°50′59. ″90. We detected the source with a count rate of (1.06 ± 0.05) × 10⁻² cts s⁻¹ in the energy range 0.2–10 keV. For the spectral analysis the 887 ± 46 pn and 363 ± 29 MOS counts were divided into 19 and 11 energy bins, respectively.

The pn and the MOS spectra were fitted simultaneously. We could not obtain acceptable fits with a single power law or blackbody, or any combination of the two (including a power law plus two blackbodies) because in all cases significant residuals were presentbelow 0.5 keV. It is very unlikely that this was due to an instrumental or calibration problem because the spectra of two nearby soft and very bright sources do not show similar features.

Excluding the energy channels below 0.5 keV, we were able to obtain a good fit with an absorbed power law with photon index Γ =3.2 ± 0.4 and N_H $=(8_{-6}^{+7})\times {10}^{20}$ cm⁻² (Fig. 1d). With N_H fixed to the DM-based value of 1.5 × 10²⁰ cm⁻² we obtained a photon index of 2.8 ± 0.1 (${\chi }_{\nu }^2=0.82$, for 18 d.o.f.). A blackbody fit was rejected (${\chi }_{\nu }^2\approx 2$ both for free and fixed N_H).

A power law fit of the pn spectrum over the whole energy range suggests the presence of two absorption lines at ~ 0.2 keV and ~0.44 keV. We therefore fitted it with a model consisting of a power law with N_H fixed to 1.5 × 10²⁰ cm⁻² and two Gaussian absorption lines. We constrained the lines to be centred at E₂ = 2 × E₁ and to have the same width σ₁ = σ₂. This model gave a good fit with E₁ = 0.222 ± 0.006 keV and ${\sigma }_1={0.015}_{-0.004}^{+0.012}$ keV (${\chi }_{\nu }^2=0.78$ for 13 d.o.f.). The best fit photon index is 2.9 ± 0.2, in good agreement with that found in the hardest part of the spectrum. The strength of the lines are poorly constrained, but they are 3σ above the 0 level (see Fig. 7). The best fit spectrum, together with its residuals and those of the single power law fit are shown in Fig. 8.

A good fit was also obtained with a power law plus blackbody with an absorption line at ~ 0.45 keV (${\chi }_{\nu }^2=0.83$ for 11 d.o.f.). However, the N_H $=({22}_{-10}^{+12})\times {10}^{20}$ cm⁻² is much higher than expected, the power law is quite steep ($\Gamma ={3.5}_{-0.5}^{+0.7}$), and the blackbody has kT = 0.045 ± 0.010 keV and a poorly constrained radius (1 < R < 70 km). This is compatible with emission from the cooling of the whole surface, but for an old pulsar such as PSR B1133+16 we would expect a lower temperature unless some reheating mechanism is operating.

To assess the statistical significance of the lines we estimated the probability of obtaining by chance a fit improvement like the observed one through Monte Carlo simulations. We simulated pn spectra of PSR B1133+16 with the same exposure time as in our observations using a model without lines, i.e. the best fit power law model with Γ = 2.36 (Table 6). We fitted each simulated spectrum with a single power law and with a power law plus two harmonically spaced lines and computed the ratio of the corresponding χ² values F = ${\chi }_{\text{PL}}^2$ / ${\chi }_{\text{PL+lines}}^2$. We found a probability of 8 × 10⁻⁴ of having a fit improvement better than that observed in the real data (i.e. F > 2.47), corresponding to a ~3.1σ significance of the lines.

For the 3D-ML spectral-timing analysis, the phases were computed using the period at the epoch of the XMM–Newton observation, P = 1.187916418694(5) s given by Szary et al. (2017). We found an average pulsed fraction of 0.29 ± 0.05, and no significant energy variation over the range 0.2–3 keV in the pn (see Fig. 9, left). A hint of pulsation was also found in the MOS, with a pulsed fraction of 0.25 ± 0.07 in the range0.2–3 keV, but it was not possible to extract the 3D spectra. The pn spectra of the total, pulsed, and unpulsed flux are plotted in Fig. 9 (right). The unpulsed spectrum cannot be fitted by a single power law (${\chi }_{\nu }^2=2.12$ for 5 d.o.f.), while it is fitted well by the same absorption features of the total model. On the contrary, the spectrum of the pulsed emission does not require the addition of the absorption lines, and a power law with $\Gamma ={2.4}_{-0.4}^{+0.5}$ satisfactorily fits the spectrum. All the results of the spectral fits are summarized in Table 6.

Fig. 6 Upper limits (3σ confidence level) on the blackbody component in PSR B0919+06, as a function of blackbody temperature. The top panel gives the limits on the emitting radius, while the lower panel gives the corresponding limits on the bolometric luminosity computed as Lbol =πR2σT4.

Fig. 7 Confidence regions (1, 2, and 3σ) of the best fit line strength vs. width of the two absorption lines in the spectrum of PSR B1133+16 (dashed lines: E1 ~ 0.22 keV, solid lines: E2 ~ 0.44 keV). The best fit values are indicated by a black plus.

Fig. 8 EPIC-pn spectrum of PSR B1133+16. The top panel shows the best fit with an absorbed power law plus two absorption lines at E1 ~ 0.22 and E2 ~ 0.44 keV. The corresponding residuals are shown in the middle panel. The bottom panel shows the residuals obtained by the best fit with a single power law in the range 0.5–10 keV.

Fig. 9 Left: pulsed fraction of PSR B1133+16 as a function of energy obtained with the 3D-ML analysis. Right: EPIC-pn total (black circles), unpulsed (blue diamonds), and pulsed (red squares) spectra of PSR B1133+16.

4 Discussion

Our ML analysis of the X-ray emission of four old rotation-powered pulsars has, in general, confirmed previous results, providing in some cases smaller uncertainties on the parameters and additional information on the spectral properties of these objects. Here, we first discuss the main results for the individual pulsars and then we put them in the broader context of the class of old rotation-powered pulsars.

In the case of PSR B0114+58, previous analysis could not distinguish between a thermal and a non-thermal spectrum (Prinz & Becker 2015), while we found that its emission is primarily thermal, being fitted well by a blackbody with kT = 0.17 ± 0.02 keV, emitting radius $R={405}_{-90}^{+110}$ m, and bolometric luminosity L_bol = (4.4 ± 3.0) × 10³⁰ erg s⁻¹. The small emitting area indicates that the thermal radiation comes from a hot spot on the stellar surface, most likely from the polar cap; in fact, the expected polar cap radius for PSR B0114+58 is $R_{\text{PC}}={((2\pi R_{*}^3)/(Pc))}^{1/2}\approx 455$ m (for a NS radius of R_* = 10 km).

We found no evidence of thermal components in the phase-averaged spectra of the other targets and derived 3σ upper limits on their bolometric luminosity of L_bol ≲ 3.2 × 10²⁸ erg s⁻¹ in PSR B0628−28 and L_bol ≲ 2.4 × 10²⁹ erg s⁻¹ in PSR B0919+06. These limits apply to temperatures in the range ~ 0.05–2 keV (see Figs. 3 and 6). However, the possible presence of a thermal component in PSR B0628−28 is suggested by the fact that its pulsed emission is softer than the unpulsed emission, and that it can be fitted by a blackbody with parameters consistent, within the errors, with the above-mentioned luminosity limit ($kT={0.10}_{-0.01}^{+0.02}$ keV, $R={150}_{-50}^{+80}$ m).

Finally, we found features below 0.5 keV with a statistical significance of ≳ 3σ in the spectrum of PSR B1133+16. This could be fitted either with a power law plus blackbody model and an absorption line at ~ 0.45 keV, or with a power law and two absorption lines at harmonically spaced energies E₁ = 0.222 keV and E₂ = 0.444 keV. The more realistic best fit parameters of the continuum component lead us to prefer the second interpretation. If the lowest energy line is an electron cyclotron feature, the implied magnetic field (B_e = 1.9 × 10¹⁰(1 + z) G, where z is the gravitational redshift) is two orders of magnitude smaller than the surface field derived from the timing parameters of the pulsar (B_d = 2.1 × 10¹² G, under the usual dipole approximation). The electrons responsible for the feature should be high in the magnetosphere, and it is unclear how they could be confined in a small region with the appropriate magnetic field value. Alternatively, if the line is attributed to protons that could be in the atmosphere close to the NS surface, the required magnetic field is B_p = 3.5 × 10¹³(1 + z) G, implying the presence of multipolar field components, as required for example in the Partially Screened Gap model (Gil et al. 2003, 2007).

4.1 Thermal X-rays from hot polar caps

In Table 7, we list the information on the non-recycled rotation-powered pulsars with evidence of thermal emission from a small region of their surface. In addition to the objects where the presence of such emission is well established (first part of the table), we have included pulsars, like PSR B0628−28, for which there is only some evidence for it (second part of the table).

The values of temperature and radius for each blackbody, with the corresponding errors, are taken from the most recent works present in literature (references in the last column of Table 7), and, when needed, we rescaled the radii using updated distance values. The bolometric luminosities of these thermal components are evaluated as L_bol = πR²σT⁴ and we included the distance uncertainties in their errors.

All the temperatures are in the narrow range 0.1–0.3 keV and do not depend significantly on the characteristic age. On the contrary, the emissionradii span nearly two orders of magnitudes. We expected the emitting region to be related to the polar cap size, so we compared R_BB with the radius of the polar cap defined by the last closed lines in a dipolar field geometry, $R_{\text{PC}}={((2\pi R_{*}^3)/(Pc))}^{1/2}$. As shown in Fig. 10, in most pulsars with 0.1 < τ < 10 Myr the two radii are consistent, considering the errors, and there is no clear correlation between R_BB/R_PC and characteristic age. Some pulsars have R_BB that is significantly smaller than R_PC, but this can be explained by geometrical effects. In fact, the radii inferred from the spectral fits correspond to the projected area of the emitting region averaged over the star rotational phase. Only for nearly aligned rotators (small angles between rotation and magnetic axis) seen pole-on does $\pi R_{\text{BB}}^2$ correspond to the real emitting area. This should occur, for example for PSR B0943+10 (Deshpande & Rankin 2001; Bilous et al. 2014); however, it has R_BB < R_PC, which might indicate that the magnetic field is not purely dipolar. The presence of multipolar components of the field can cause a reduction of the polar cap area by a factor B_s∕B_d, where B_s and B_d are the actual magnetic field at the star surface and that of the dipole, respectively (Gil & Sendyk 2000). PSR B0656+14, on the other hand, has a blackbody radius that is nearly ten times larger than R_PC, as already noticed by De Luca et al. (2005), suggesting the need for a different explanation of the hot thermal component in this pulsar. One possibility is that anisotropic thermal conduction in the crust plays a role in causing temperature gradients on the surface, andthe oversimplified modelling with just two blackbodies at different temperatures gives unrealistic parameters.

It is also interesting to compare the thermal polar cap emission of old (non-recycled) pulsars with that observed in millisecond pulsars. Thermal X-rays have been detected in several millisecond pulsars (Bogdanov et al. 2006; Zavlin 2006; Spiewak et al. 2016; Bhattacharya et al. 2017). For two of them (PSR J0030+0451 Bogdanov & Grindlay 2009; PSR J0437−4715 Bogdanov 2013) at least two thermal components are needed, which are interpreted as emission from the polar cap having a non-uniform temperature distribution with a hot core and a cooler rim. We have included millisecond pulsars in Fig. 10 where we used the emission radii derived from blackbody fits in order to make a consistent comparison with the values of non-recycled pulsars. We note that fits with a non-magnetized hydrogen atmosphere model, which are appropriate for millisecond pulsars, would result in emitting radii larger by a factor ~4, thus making most of them consistent with the polar cap radii in the centred dipole approximation.

Thermal emission from small hot spots has been explored in detail by Harding & Muslimov (1998, 2001, 2002) in the context of polar cap models. In their Space Charge Limited Flow model, primary particles accelerated above the polar cap produce a photon-pairs cascade in an avalanche process. Backward particles screen the acceleration voltage drop and heat the polar cap. The pair cascade can be initiated either by curvature radiation (CR) or by inverse-Compton scattering (ICS) of thermal X-rays by primary electrons. The production of CR photons requires a much higher energy of the primary particles, and therefore the pair production occurs at higher altitude and the pair heating luminosity is expected to be much higher in the CR rather than in the ICS scenario. These processes can occur only above the corresponding “death lines” in the Ṗ −P diagram. Pulsars in the region below the CR pairs death line, but above that of ICS, can produce pairs only from ICS photons.

Harding & Muslimov (2002) give analytic expressions for the expected surface X-ray luminosity, depending on the pair fronts creation mechanism. In the case of non-recycled pulsars, L_bol ≈ 10³¹ erg s⁻¹ P^−1∕2 for CR and L_bol ≈ 2.0 × 10²⁷ erg s⁻¹ P^−3∕2 for ICS. These relations, which should be considered as upper limits in the case of incomplete screening, are compared with the observed thermal luminosities of the pulsars in Fig. 11. Although there is a general agreement between the data and the model predictions, we note that among the three pulsars for which only the ICS mechanism should operate (indicated in Fig. 11 by a square), two have luminosities above the maximum predicted in this scenario. One of them is PSR B0628−28, for which the possible thermal origin of the pulsed flux is admittedly not very compelling. On the other hand, the presence of a thermal component in PSR B0943+10 is well established (Mereghetti et al. 2016) and possibly requires a different explanation.

Fig. 10 Ratio of the radius of emitting area inferred from the blackbody fit to the polar cap radius as a function of the pulsar characteristic age. Millisecond pulsars are indicated with crosses. In the case of multiple thermal components (PSR J0030+0451 two components, PSR J0437−4715 three components), all the radii are plotted.

Fig. 11 Bolometric luminosity of hot polar cap emission versus rotation period. The lines represent the predicted thermal luminosities in case of curvature radiation (red dashed line) and inverse-Compton scattering (blue dotted line) pair fronts. The theoretical lines are from Harding & Muslimov (2002) in the case of non-recycled pulsars. Pulsars marked with a square are expected to produce pairs only from ICS photons. The error bars on the bolometric luminosity are computed taking into account both the uncertainties of the blackbody fit and those on the distances. See Table 7 for details.

5 Conclusions

We applied the powerful ML technique to the spectral X-ray analysis of four dim rotation-powered old pulsars and found evidence for absorption lines at 0.22 and 0.44 keV in the spectrum of PSR B1133+16. If these lines are interpreted as cyclotron features due to protons they imply a magnetic field of a few 10¹³ G, higher thanthe dipole field inferred from the timing parameters of the pulsar (B_d= 2 × 10¹² G).

The presence of thermal emission from a hot polar cap could be established only for PSR B0114+58. Data of better statistical quality are required to ascertain if thermal components are really needed, also in the three other pulsars considered here: PSR B0628−28, PSR B0919+06, and PSR B1133+16. Several old pulsars show the presence of thermal components in their X-ray spectra. By using the most updated spectral results and distance estimates, we found that the sizes of the emitting hot spots in these pulsars are generally consistent with the expected dimensions of the polar caps once projection effects are taken into account.

Acknowledgements

This work has been partially supported through PRIN INAF 2014.

References

All Tables

All Figures

	Fig. 1 Confidence regions (1, 2, and 3σ) of the best fit parameters for the spectra of (a) PSR B0628−28, (b) PSR B0919+06, (c) PSR B0114+58, and (d) PSR B1133+16. The best fit values are indicated by a plus and the DM-based value of NH by a dashed line.
In the text

	Fig. 2 EPIC-pn (blue diamonds) and -MOS (red squares) X-ray spectra of PSR B0628−28. The best fit power law is shown in the top panel, and the corresponding residuals in the middle panel. The bottom panel shows the residuals obtained by fitting the spectra with an absorbed blackbody.
In the text

	Fig. 3 Upper limits (3σ confidence level) on the blackbody component in PSR B0628−28 as a function of blackbody temperature. The top panel gives the limits on the emitting radius, while the lower panel gives the corresponding limits on the bolometric luminosity computed as Lbol =πR2σT4.
In the text

	Fig. 4 Left: pulsed fraction of PSR B0628−28 as a function of energy obtained with the 3D-ML analysis. Right: EPIC-pn total (black circles), unpulsed (blue diamonds), and pulsed (red squares) spectra of PSR B0628−28.
In the text

	Fig. 5 EPIC-pn (blue diamonds) and -MOS (red squares) X-ray spectra of PSR B0919+06. The best fit power law is shown in the top panel, and the corresponding residuals in the middle panel. The bottom panel shows the residuals obtained by fitting the spectra with an absorbed blackbody.
In the text

	Fig. 6 Upper limits (3σ confidence level) on the blackbody component in PSR B0919+06, as a function of blackbody temperature. The top panel gives the limits on the emitting radius, while the lower panel gives the corresponding limits on the bolometric luminosity computed as Lbol =πR2σT4.
In the text

	Fig. 7 Confidence regions (1, 2, and 3σ) of the best fit line strength vs. width of the two absorption lines in the spectrum of PSR B1133+16 (dashed lines: E1 ~ 0.22 keV, solid lines: E2 ~ 0.44 keV). The best fit values are indicated by a black plus.
In the text

	Fig. 8 EPIC-pn spectrum of PSR B1133+16. The top panel shows the best fit with an absorbed power law plus two absorption lines at E1 ~ 0.22 and E2 ~ 0.44 keV. The corresponding residuals are shown in the middle panel. The bottom panel shows the residuals obtained by the best fit with a single power law in the range 0.5–10 keV.
In the text

	Fig. 9 Left: pulsed fraction of PSR B1133+16 as a function of energy obtained with the 3D-ML analysis. Right: EPIC-pn total (black circles), unpulsed (blue diamonds), and pulsed (red squares) spectra of PSR B1133+16.
In the text

	Fig. 10 Ratio of the radius of emitting area inferred from the blackbody fit to the polar cap radius as a function of the pulsar characteristic age. Millisecond pulsars are indicated with crosses. In the case of multiple thermal components (PSR J0030+0451 two components, PSR J0437−4715 three components), all the radii are plotted.
In the text

Fig. 11 Bolometric luminosity of hot polar cap emission versus rotation period. The lines represent the predicted thermal luminosities in case of curvature radiation (red dashed line) and inverse-Compton scattering (blue dotted line) pair fronts. The theoretical lines are from Harding & Muslimov (2002) in the case of non-recycled pulsars. Pulsars marked with a square are expected to produce pairs only from ICS photons. The error bars on the bolometric luminosity are computed taking into account both the uncertainties of the blackbody fit and those on the distances. See Table 7 for details.

In the text