XMM-Newton observations of PSR J0726-2612, a radio-loud XDINS

TL;DR

This study uses XMM-Newton data to analyze PSR J0726-2612, revealing thermal X-ray emission, a broad absorption line, and a pulse profile suggesting magnetic beaming, thereby strengthening its similarity to X-ray Dim Isolated Neutron Stars (XDINSs).

Contribution

First detailed X-ray spectral and timing analysis of PSR J0726-2612, confirming its thermal emission and magnetic beaming, and comparing its geometry with XDINSs.

Findings

Thermal X-ray emission fitted by two blackbodies.

Presence of a broad absorption line at 0.39 keV.

Pulse profile indicates magnetic beaming effects.

Abstract

We present the results of an XMM-Newton observation of the slowly rotating ($P = 3.4$ s), highly magnetized ($B \approx 3 \times 10^{13}$ G) radio pulsar PSR J0726-2612. A previous X-ray observation with the Chandra satellite showed that some of the properties of PSR J0726-2612 are similar to those of the X-ray Dim Isolated Neutron Stars (XDINSs), a small class of nearby slow pulsars characterized by purely thermal X-ray spectra and undetected in the radio band. We confirm the thermal nature of the X-ray emission of PSR J0726-2612, which can be fit by the sum of two blackbodies with temperatures $kT_1 = 0.074_{-0.011}^{+0.006}$ keV and $kT_2 = 0.14_{-0.02}^{+0.04}$ keV and emitting radii $R_1 = 10.4_{-2.8}^{+10.8}$ km and $R_2 = 0.5_{-0.3}^{+0.9}$ km, respectively (assuming a distance of 1 kpc). A broad absorption line modeled with a Gaussian profile centred at $0.39_{-0.03}^{+0.02}$ keV is required in the fit. The pulse profile of PSR J0726-2612 is characterized by two peaks with similar intensity separated by two unequal minima, a shape and pulsed fraction that cannot be reproduced without invoking magnetic beaming of the X-ray emission. The presence of a single radio pulse suggests that in PSR J0726-2612 the angles that the dipole axis and the line of sight make with the rotation axis, $\xi$ and $\chi$ respectively, are similar. This geometry differs from that of the two radio-silent XDINSs with a double peaked pulse profile similar to that of PSR J0726-2612, for which $\xi \sim 90^\circ$ and $\chi \sim 45^\circ$ have been recently estimated. These results strengthen the similarity between PSR J0726-2612 and the XDINSs and support the possibility that the lack of radio emission from the latter might simply be due to an unfavourable viewing geometry.

Tables

Table 1. Table 1: Observed and derived parameters for PSR J0726 − - 2612

R.A. (J2000.0)	${07}^h{26}^m{08}^s.12\left(4\right)$
Dec. (J2000.0)	$-{26}^{\circ }{12}^{\prime }{38}^{\prime \prime }.1\left(8\right)$
Period $P$ (s)	$3.4423084877\left(4\right)$
Period derivative $\dot{P}$ (s s-1)	$2.9311\left(4\right)\times {10}^{-13}$
Epoch (MJD)	$52,950$
Characteristic age ${\tau }_c$ (years)	$1.86\times {10}^5$
Surface dipolar magnetic field $B_s$ (G)	$3.2\times {10}^{13}$
Rotational energy loss rate $\dot{E}$ (erg s-1)	$2.8\times {10}^{32}$
Dispersion measure DM (cm-3 pc)	$69.4\left(4\right)$

Table 2. Table 2: Exposure Times and Source Counts for PSR J0726 − - 2612 in the three EPIC cameras

Data	EPIC camera	Exposure time	Source Counts
		ks	$0.15-1.5$ keV
Phase-averaged	pn	37.8	$18,938\pm 140$
	MOS1	64.0	$4,499\pm 69$
	MOS2	70.4	$5,212\pm 74$
Min 1	pn	9.4	$3,823\pm 63$
Max 1	pn	9.4	$5,576\pm 76$
Min 2	pn	9.4	$4,088\pm 65$
Max 2	pn	9.4	$5,447\pm 75$

Table 3. Table 3: Results for the phase-averaged and phase-resolved spectra of PSR J0726 − - 2612

Model	$N_H$ a	$k{T}_1$	$R_1^b$	$k{T}_2$	$R_2^b$	$E$	$\sigma$	strength${}^c$	$F_{\mathrm{unabs}}^{0.1-2}$	${\chi }_{\nu }^2$/dof	nhp
	${10}^{20}$ cm-2	keV	km	keV	km	keV	keV	keV	erg s-1 cm-2
Phase-averaged spectra:
BB	$4.1\pm 0.2$	$0.0896\left(6\right)$	$4.90\pm 0.15$	…	…	…	…	…	${1.60}_{-0.05}^{+0.06}$	1.37/213	$3\times {10}^{-4}$
2BB	$4.3\pm 0.2$	$0.0888\left(7\right)$	$5.1\pm 0.2$	$>0.33$		…	…	…	${1.65}_{-0.06}^{+0.07}$	1.32/211	$1\times {10}^{-3}$
GBB	$2.8\pm 0.3$	$0.11\pm 0.01$	${2.9}_{-0.4}^{+0.5}$	…	…	$1.09\pm 0.09$	$0.28\pm 0.08$	${1.0}_{-0.6}^{+1.1}$	${1.37}_{-0.09}^{+0.35}$	1.12/210	$0.11$
G2BB	${5.3}_{-0.8}^{+1.2}$	${0.074}_{-0.011}^{+0.006}$	${10.4}_{-2.8}^{+10.8}$	${0.14}_{-0.02}^{+0.04}$	${0.5}_{-0.3}^{+0.9}$	${0.39}_{-0.03}^{+0.02}$	${0.08}_{-0.02}^{+0.03}$	${0.12}_{-0.05}^{+0.13}$	${3.30}_{-0.85}^{+3.85}$	1.00/208	$0.47$
GNSA${}^e$	${6.9}_{-1.1}^{+0.8}$	$0.029\left(1\right)$	${14.3}^d$	…	…	${0.37}_{-0.03}^{+0.02}$	${0.09}_{-0.01}^{+0.02}$	${0.17}_{-0.04}^{+0.07}$	$9.0\pm 1.1$	1.03/210	$0.36$
GNSMAXG${}^f$	${5.9}_{-4.2}^{+3.4}$	$0.028\left(1\right)$	${14.3}^d$	…	…	$0.28\pm 0.09$	${0.14}_{-0.04}^{+0.06}$	${0.62}_{-0.33}^{+1.26}$	${19.3}_{-7.0}^{+9.7}$	1.02/210	$0.40$
G2BB phase-resolved:
Maxima 1	${5.3}^d$	${0.074}^d$	${10.4}^d$	$0.110\left(7\right)$	${1.55}_{-0.30}^{+0.40}$	$0.39\pm 0.01$	$0.07\pm 0.01$	${0.11}_{-0.01}^{+0.02}$	$3.6\pm 0.2$	1.00/80	$0.47$
Maxima 2	${5.3}^d$	${0.074}^d$	${10.4}^d$	$0.111\left(9\right)$	${1.40}_{-0.35}^{+0.45}$	$0.39\pm 0.01$	$0.06\pm 0.01$	$0.10\pm 0.01$	$3.6\pm 0.2$	0.95/78	$0.61$
Minima 1	${5.3}^d$	${0.074}^d$	${10.4}^d$	${0.17}_{-0.04}^{+0.07}$	${0.20}_{-0.13}^{+0.31}$	${0.40}_{-0.02}^{+0.01}$	${0.14}_{-0.02}^{+0.04}$	${0.23}_{-0.03}^{+0.05}$	${2.90}_{-0.25}^{+0.45}$	1.14/49	$0.23$
Minima 2	${5.3}^d$	${0.074}^d$	${10.4}^d$	${0.29}_{-0.09}^{+0.72}$	${0.06}_{-0.03}^{+0.05}$	${0.39}_{-0.03}^{+0.02}$	${0.13}_{-0.02}^{+0.03}$	${0.20}_{-0.02}^{+0.04}$	${3.1}_{-0.2}^{+0.3}$	1.28/52	$0.08$

Table 4. Table 4: Comparison between PSR J0726 − - 2612 and the XDINSs

Source	$P$	$\dot{P}$	$B_p$	$E_{\mathrm{cyc}}$	$B_{\mathrm{cyc}}$	$L_X/\dot{E}$	Pulse	PF	Refs.
RX	s	${10}^{-14}$ s s-1	${10}^{13}$ G	eV	${10}^{13}$ G			%
J0420.0$-$5022	3.45	2.76	2.0	…	…	$0.31-0.38$	single	13	(1)
J0720.4$-$3125	16.78	18.6	11.3	${254}_{-30}^{+25}$	${3.4}_{-0.4}^{+0.3}$	$99-157$	double	11	(2)
J0806.4$-$4123	11.37	5.6	5.1	${241}_{-12}^{+11}$	${3.2}_{-0.2}^{+0.2}$	$10.6-16.7$	single	6	(1)
J1308.6$+$2127 (RBS 1223)	10.31	11.2	6.9	${390}_{-6}^{+6}$	${5.16}_{-0.08}^{+0.08}$	$31.5-39.6$	double	18	(3)
J1605.3$+$3249 (RBS 1556)	…	…	…	${353}_{-48}^{+19}$	${4.7}_{-0.6}^{+0.3}$	…	…	1.4	(4)
J1856.5$-$3754	7.06	2.98	2.9	…	…	$9.6-15.2$	single	1.2	(5)
J2143.0$+$0654 (RBS 1774)	9.43	4.1	4.0	${326}_{-79}^{+56}$	${4.3}_{-1.0}^{+0.7}$	$33.2-41.8$	single	4	(6-7)
PSR J0726$-$2612	3.44	29.3	6.4	${390}_{-20}^{+10}$	${5.2}_{-0.3}^{+0.1}$	$1.1-3.0$	double	30	(8)

Full text

††offprints: m.rigoselli@campus.unimib.it, michela.rigoselli@inaf.it

1 Introduction

Observations with the ROSAT satellite in the mid-1990s led to the discovery of a small group of isolated neutron stars characterized by thermal emission at soft X-rays, now known as XDINSs (X-ray Dim Isolated Neutron Stars, see Haberl 2007; Turolla 2009 for reviews). XDINSs have spin periods in the range $P\sim 3-17$ s and period derivatives of a few 10*-14* s s*-1*, which result in characteristic ages $\tau_{c}=P/2\dot{P}\sim 1-4$ Myr. With the usual assumption that the spin-down is due to magnetic dipole braking, these timing parameters imply magnetic fields of the order of a few $10^{13}$ G.

The XDINSs are at distances of only a few hundreds parsecs and for two of them the parallax of the optical counterpart has been measured (Walter et al. 2010; Tetzlaff et al. 2011). The XDINSs have X-ray luminosities of $10^{31}-10^{32}$ erg s*-1*, higher than their spin-down power. Their X-ray spectra are very soft, with blackbody temperatures of $kT\sim 45-110$ eV, often showing the presence of broad absorption lines. If these lines are interpreted as proton cyclotron features or atomic transitions (see, e.g., Kaplan 2008), the magnetic fields estimated from their energies are of the same order of those derived from the spin-down rate assuming magnetic dipole braking. The X-ray emission of XDINSs, consisting only of thermal components, is believed to come directly from the star surface and, given the small distance of these sources, it is little affected by photoelectric absorption in the interstellar medium. The XDINS discovery raised some excitement since they appeared as optimal targets to test neutron star surface emission models without being affected by the presence of non-thermal emission. However, the ultimate goal of constraining the star radius and hence the equation of state with these studies is still hampered by our poor knowledge of the neutron star surface layers composition and magnetization.

The attempt to explain the different manifestations of neutron stars (e.g. Mereghetti 2011) in the context of a unified evolutionary picture is one of the current challenges in the study of neutron stars (Kaspi 2010; Igoshev et al. 2014). In the $P-\dot{P}$ diagram, shown in Fig. 1, XDINSs are located in the region below that occupied by the magnetars, a group of isolated neutron stars powered mainly by magnetic energy (see, e.g. Mereghetti et al. 2015; Turolla et al. 2015; Kaspi & Beloborodov 2017). This has led to the suggestion that the XDINSs might be the descendent of magnetars (Heyl & Kulkarni 1998; Colpi et al. 2000). The strong internal field of magnetars ($B\gtrsim 10^{15}$ G) significantly affects their thermal evolution (Viganò et al. 2013), resulting in luminosities higher than those predicted for normal pulsars of similar age.

A distinctive property of the XDINSs is that they are not detected in the radio band111The possible detection of pulsed emission from two XDINSs at very low frequencies (Malofeev et al. 2005, 2006) is, so far, unconfirmed. (Kondratiev et al. 2009). The reason for the lack of radio emission is still uncertain. One possibility is that this is due to their old age and long spin period (Baring & Harding 1998, 2001). However, a few radio pulsars with periods $\gtrsim\!10$ s have been recently discovered: PSR J0250$+$5854 with $P=23.5$ s (Tan et al. 2018), and a second one with $P=12.1$ s (Morello et al. 2019 in prep). Another explanation might be related to the geometrical configuration of their magnetosphere, that, especially if they are old magnetars, might be strongly non-dipolar (Turolla et al. 2015). Finally, it cannot be excluded that (at least some of) the XDINSs are simply ordinary radio pulsars with radio beams unfavorably aligned with respect to the Earth. In this respect, it is interesting to investigate radio-loud pulsars with X-ray properties and/or timing parameters similar to those of the XDINSs, such as the long period (greater than a few seconds) and high B ($>10^{13}$ G) pulsars.

Among these, here we focus on PSR J0726$-$2612, a radio pulsar with spin period $P=3.44$ s and characteristic age of 200 kyr that was discovered in the Parkes High-Latitude Survey (Burgay et al. 2006). Its timing parameters (Table 1) are in the range of those of the XDINSs. The similarity with the XDINSs was reinforced by X-ray observations with the Chandra satellite (Speagle et al. 2011), that revealed a soft thermal spectrum with blackbody temperature $kT\approx 87$ eV, and pulsations with a sinusoidal, double-peaked profile. The distance of PSR J0726$-$2612 is unknown. Its dispersion measure DM $=69.4\pm 0.4$ cm*-3* pc (Burgay et al. 2006) implies a distance $d=2.9$ kpc, assuming the Galactic electrons distribution of Yao et al. (2017). However, there are a few facts suggesting that this is probably an overestimate. Such a large value would give a distance of 230 pc from the Galactic plane, implying, if PSR J0726$-$2612 was born close to the plane and its true age is similar to $\tau_{c}$, a velocity of the order of a thousand km s*-1*. This value is not impossible, but it would be at the far end of the pulsar velocity distribution (Hobbs et al. 2005). More importantly, for such a large $d$, one would expect an X-ray absorption corresponding to a sizeable fraction of the total Galactic H I column density, that in this direction is $\sim 5\times 10^{21}$ cm*-1* (Kalberla et al. 2005), while the observed value is a factor 10 smaller. Finally, the line of sight toward PSR J0726$-$2612 crosses the Gould belt, that is not included in the electron distribution model of Yao et al. (2017). This could explain the large distance inferred from the DM. This local structure ($d\sim 200-400$ pc) comprises several OB associations that have been proposed as the birthplace of the XDINSs (Popov et al. 2003, 2005). Speagle et al. (2011) suggested that also PSR J0726$-$2612 could be associated with the Gould belt and hence closer than $\sim 1$ kpc.

Here we report the results of XMM–Newton observations which show other similarities between PSR J0726$-$2612 and the XDINSs. In the following we will scale all the distance-dependent quantities to $d_{\rm kpc}=1$ kpc and adopt representative values of mass and radius of $1.2\leavevmode\nobreak\ M_{\odot}$ and 12 km, respectively.

2 Observations and data reduction

PSR J0726$-$2612 was observed with the European Photon Imaging Cameras (EPIC) instrument on board XMM–Newton with a single pointing lasting 108 ks on 2013 April 8. The three cameras of EPIC ($0.1-12$ keV), the pn (Strüder et al. 2001) and the two MOS (Turner et al. 2001), were operated in Full Frame mode with the thin optical filter. While the pn time resolution (73.4 ms) is adequate to reveal the pulsations of the source, this is impossible for the MOS given its resolution time of 2.6 s.

The data reduction was performed using the epproc and emproc pipelines of version 15 of the Science Analysis System (SAS)222https://www.cosmos.esa.int/web/xmm-newton/sas. We selected single- and multiple-pixel events (PATTERN $\leq 4$ and PATTERN $\leq 12$) for both the pn and MOS. We then removed time intervals of high background using the SAS program espfilt with standard parameters. The source was detected by EPIC at coordinates R.A. = ${07}^{\rm h}{26}^{\rm m}{08}^{\rm s}{.1}$, Dec. = ${-26}^{\circ}{12}^{\prime}{38}^{\prime\prime}{}$, fully consistent with the radio position (Table 1). The source events were selected from a circle of radius $40^{\prime\prime}$ centred at the radio position, while the background was extracted from a nearby circular region of radius $60^{\prime\prime}$. The resulting net exposure times and source events are listed in Table 2. At the corresponding count rates pile-up effects are not relevant.

3 Results

3.1 Timing analysis

PSR J0726$-$2612 is barely detected above 1.5 keV, therefore we limited our timing analysis to the energy band $0.15-1.5$ keV. The times of arrival were converted to the barycenter of the Solar System with the task barycen. An epoch folding search of the EPIC-pn data gave a best period $P=3.442396(1)$ s, that is consistent within $0.7\sigma$ with the value expected at the XMM–Newton observation epoch ($56,390$ MJD) using the ATNF ephemeris reported in Table 1. The background-subtracted light curve, in the energy band $0.15-1.5$ keV is shown in Fig. 2. The position of the radio pulse is indicated, with its 1$\sigma$ uncertainty, as a vertical red line.

The EPIC-pn pulse profile shows two peaks with the same intensity (net count rate of max${}_{1}=0.62\pm 0.02$ cts s*-1* and max${}_{2}=0.64\pm 0.02$ cts s*-1*), separated by about 0.5 cycles. The two minima of the pulse profile are instead significantly different: min${}_{1}=0.34\pm 0.01$ cts s*-1* and min${}_{2}=0.39\pm 0.01$ cts s*-1*. The pulse profile is symmetric in phase with respect to any of the two minima, but a fit with a constant plus a sine function at half of the spin period is not acceptable ($\chi^{2}_{\nu}=2.7$ for 17 dof). The pulsed fraction333Defined as (max(CR)-min(CR))/(max(CR)+min(CR)), where CR is the background-subtracted count rate. is $30\pm 2\%$.

Fig. 3 shows that the soft ($0.15-0.4$ keV) and hard energy ranges ($0.4-1.5$ keV) have slightly different pulsed fractions: $26\pm 3\%$ and $37\pm 3\%$, respectively. Moreover, the positions of the first minimum and of the second maximum are shifted of about 1 bin between the two energy ranges, but the symmetry around the minima is preserved in both bands. Fits with a constant plus sine function give $\chi^{2}_{\nu}=1.6$ and $\chi^{2}_{\nu}=3.5$ for the soft and hard profile, respectively. The hardness ratio444Defined as (hard(CR)-soft(CR))/(hard(CR)+soft(CR)), where the soft energy range is $0.15-0.4$ keV, the hard one $0.4-1.5$ keV., shown in the lower panel of the same figure, clearly indicates the presence of phase-dependent spectral variations: the source is softer during the minima and harder during the maxima.

3.2 Spectral analysis

The spectral analysis was performed using XSPEC (ver. 12.8.2). The spectra were rebinned using the GRPPHA tool with a minimum of 50 counts per bin. The spectra of the three cameras were fitted simultaneously, including a renormalization factor to account for possible cross-calibration uncertainties. Errors on the spectral parameters are at $1\sigma$ confidence level.

We used the photoelectric absorption model tbabs, with cross sections and abundances from Wilms et al. (2000). Both a single power law and a blackbody did not provide acceptable fits, giving $\chi^{2}_{\nu}\approx 6$ and $\chi^{2}_{\nu}=1.37$ for 213 dof (Null Hypothesis Probability, nhp, of $3\times 10^{-4}$), respectively. We then attempted a fit with magnetized hydrogen atmosphere models (nsa and nsmaxg in XSPEC, Pavlov et al. 1995; Ho et al. 2008; Ho 2014). However, none of the two sets of available models (the first with a single surface $B$ and $T_{\rm eff}$, the second with $B$ and $T_{\rm eff}$ varying across the surface according to the magnetic dipole model) gave an acceptable fit ($\chi^{2}_{\nu}>2.2$ for 213 dof). In conclusion, we could not find a good fit with single component models.

Also modelling the spectra with a blackbody plus power law or with the sum of two blackbodies was unsatisfactory. In the first case we obtained a negative photon index for the power law, while in the second case, the second thermal component had a negligible flux, and did not improve the quality of the fit with respect to that of a single blackbody ($\chi^{2}_{\nu}=1.32$ for 211 dof, nhp $=10^{-3}$).

A real improvement in the fit was obtained by adding to the blackbody a broad absorption line modelled with a Gaussian (GBB) centered at $E=1.09\pm 0.09$ keV and width $\sigma=0.28\pm 0.08$ keV ($\chi^{2}_{\nu}=1.12$ for 210 dof). Following the recent results of Yoneyama et al. (2019), we explored the possibility to adopt a two blackbody component model plus a Gaussian line in absorption (G2BB). With this model we found a good fit with the line placed at $E=0.39_{-0.03}^{+0.02}$ keV and with a broadening of $\sigma=0.08_{-0.02}^{+0.03}$ keV ($\chi^{2}_{\nu}=1.00$ for 208 dof). The addition of the line yields an improvement of the $\chi^{2}$ of $F=\chi^{2}_{\rm 2BB}/\chi^{2}_{\rm G2BB}=1.32$. To assess the statistical significance of the line, we estimated through Monte Carlo simulations the probability of obtaining by chance an equal (or better) fit improvement: we estimate a probability of $\sim 10^{-5}$ of having $F\geq 1.32$, corresponding to a $\sim 4.4\sigma$ significance of the line. The cold blackbody ($kT_{1}\approx 0.074$ keV) has an emitting radius $R_{1}=10.4_{-2.8}^{+10.8}$ $d_{\rm kpc}$ km, compatible with emission from the whole neutron star, while the hot blackbody has $kT_{2}\approx 0.14$ keV and $R_{2}=0.5_{-0.3}^{+0.9}$ $d_{\rm kpc}$ km.

A good fit was also found with the magnetized atmosphere models with a dipole distribution of the surface magnetic field ($B=10^{13}$ G at the poles) plus a Gaussian line in absorption. With the nsa model, we found an effective temperature $T_{\rm eff}=0.40\pm 0.08$ MK (corresponding to an observed temperature $kT=0.029\pm 0.001$ keV), $d=121_{-12}^{+13}$ pc and $E=0.37_{-0.03}^{+0.02}$ keV, $\sigma=0.09_{-0.01}^{+0.02}$ keV for the Gaussian line ($\chi^{2}_{\nu}=1.03$ for 210 dof). With the nsmaxg model, for an impact parameter (that is the angle between the line of sight and the dipole axis) $\eta=90^{\circ}$, the model parameters are $T_{\rm eff}=0.39\pm 0.02$ MK ($kT=0.028\pm 0.001$ keV), $d=63_{-17}^{+26}$ pc and $E=0.28\pm 0.09$ keV, $\sigma=0.14_{-0.04}^{+0.06}$ keV for the Gaussian line ($\chi^{2}_{\nu}=1.02$ for 210 dof). Using instead the same model with $\eta=0^{\circ}$, the fit was not acceptable ($\chi^{2}_{\nu}=2.38$ for 210 dof).

The spectral results are summarized in Table 3, while in Fig. 4 the best blackbody fits are shown.

The light curves and hardness ratio shown in Fig. 3 indicate that a spectral variation occurs as a function of the rotation phase. Therefore, we extracted the EPIC-pn spectra of the phase intervals corresponding to the two minima and the two maxima of the pulse profile, as shown in Fig. 2 (the number of source events in each spectrum is listed in Table 2). In order to illustrate the spectral variations, we fitted the spectra with the G2BB model, fixing all of the parameters at the best fit values of the phase-averaged spectrum, except for an overall normalization. The residuals, shown in the two lower panels of Fig. 5, indicate that the spectra of the two maxima are similar and significantly harder then those of the minima. Their normalization factors with respect to the phase-averaged spectrum, $N_{\rm max1}=1.16\pm 0.02$ and $N_{\rm max2}=1.14\pm 0.02$, are consistent with the same value, while those of the two minima are different ($N_{\rm min1}=0.85\pm 0.01$ and $N_{\rm min2}=0.80\pm 0.01$).

We then fitted the four spectra separately, keeping fixed only the interstellar absorption and the parameters of the cold blackbody, because we do not expect them to vary during a stellar rotation. The results are given in Table 3. The absorption line is at the same energy in the four spectra, but it has different widths and normalizations. The hot blackbody temperature is lower ($kT\approx 0.11$ keV) and its emission radius is larger ($R\approx 1.5$ $d_{\rm kpc}$ km) at the two maximum phases than at the first minimum ($kT\approx 0.17$ keV and $R\approx 0.20$ $d_{\rm kpc}$ km), while these parameters are poorly constrained at the second minimum. We also tried other fits allowing more parameters to vary, but the results were inconclusive due to the strong parameter degeneracy.

4 Discussion

Our XMM–Newton results for PSR J0726$-$2612 are consistent with those previously obtained with Chandra (Speagle et al. 2011), but, thanks to a significant detection with good statistics over a broader energy range, they provide more information on the spectrum and pulse profile of this pulsar.

We found that the spectrum of PSR J0726$-$2612 is more complex than the single blackbody that was adequate to fit the Chandra data. The single blackbody fit requires the addition of a broad absorption line at $E\approx 1.09$ keV. A better fit was obtained with two blackbody components, but also in this case a line at $E\approx 0.39$ keV is required. The colder blackbody component has an emitting area consistent with a large fraction of the star surface ($R_{1}=10.4_{-2.8}^{+10.8}$ $d_{\rm kpc}$ km), while the hotter one can be attributed to a small hot spot ($R_{2}=0.5_{-0.3}^{+0.9}$ $d_{\rm kpc}$ km), likely located at the magnetic pole.

Our results confirmed that the interstellar absorption is about a factor of ten smaller than the value ($N_{\rm H}$ $=2.1\times 10^{21}$ cm*-2*) inferred from the dispersion measure and the usual assumption of a $10\%$ ionization of the interstellar medium (He et al. 2013). This might be due to the line of sight crossing the Gould belt.

An equally good fit was obtained with a magnetized hydrogen atmosphere covering all the star surface, but also in this case the presence of an absorption line at $E\approx 0.37$ keV (nsa model) or $E\approx 0.28$ keV (nsmaxg model) is required. We note that the constant (polar) value of the magnetic field in the nsa (nsmaxg) model is fixed in the fits at $B=10^{13}$ G, and that the nsa model assumes a uniform distribution of the temperature. The nsmaxg model is more realistic, but it assumes that the dipole axis is orthogonal to the line of sight, that does not necessarily apply to the case of PSR J0726$-$2612. Moreover, the inferred distance of $\approx 63$ pc seems unrealistically small.

The absorption lines we found in the spectra can be interpreted as proton cyclotron features at $E_{\rm cyc}=0.063\leavevmode\nobreak\ B_{13}\times(1+z)$ keV, where $z$ is the gravitational redshift and $B_{13}$ the magnetic field in units of $10^{13}$ G. In the case of G2BB model, for $E_{\rm cyc}=0.39$ keV and $z\approx 0.2$, we get $B\approx 5\times 10^{13}$ G, in good agreement with the dipole magnetic field evaluated at the poles ($B_{p}\approx 6\times 10^{13}$ G). However, we caution that other explanations cannot be ruled out, including the possibility that the lines are simply an artefact resulting from an oversimplified modeling of the continuum emission. In fact, Viganò et al. (2014) showed that non-homogeneus temperature distributions on a neutron star surface can, in some cases, lead to the appearance of broad features when the spectra are fitted with simple blackbody models.

Contrary to the previous Chandra results, we also found that the double-peaked pulse profile of PSR J0726$-$2612 is not well described by a sinusoid, owing to the significant difference in the flux of the two minima. Remarkably, the pulse profile is symmetric for phase reflection around any of the two minima. Within the limits due to their lower statistics, these properties seem to hold also for the profiles in the soft and hard X-ray bands. The pulse profiles are moderately energy-dependent, with evidence for a harder emission in correspondence of the two peaks.

Although a detailed modeling of the light curves of PSR J0726$-$2612 is beyond the scope of the present work, we explored whether a simple model based on blackbody emission components with parameters consistent with the spectral results could reproduce the pulse profile. We assumed that the hotter blackbody comes from two antipodal magnetic polar caps with opening angle $\theta_{c1}=3^{\circ}$, while the colder one from two annuli extending between $\theta_{c1}$ and $\theta_{c2}=36^{\circ}$. The temperatures of the emitting regions were set to the values derived from the spectral analysis (model G2BB, $kT_{1}=0.074$ keV, $kT_{2}=0.14$ keV) and the angular apertures were chosen in such a way to reproduce the emitting radii derived from the fit for a NS radius of 12 km. We also added interstellar absorption and a Gaussian absorption line, with parameters fixed to those of the phase averaged spectrum. Synthetic light curves were computed using the method by Turolla & Nobili (2013) and account for general-relativistic effects. We convolved the obtained light curves with the EPIC-pn instrumental response and we evaluated the pulsed fraction in the energy range $0.15-1.5$ keV. The results depend on the angles $\chi$ and $\xi$ that the rotation axis makes with the line of sight and the magnetic axis, respectively. As shown in Fig. 6, this simple model is unable to yield the observed pulsed fraction even for the most favourable geometry (PF $\approx 21\%$ for $\xi\approx\chi\gtrsim 35^{\circ}$). This is also true if only two antipodal point-like polar caps are considered, which is the configuration yielding the maximum pulsed fraction using isotropic emission (see e.g. Turolla & Nobili 2013). Another problem is that, owing to the intrinsic symmetry of the model, the resulting light curves cannot exhibit different minima, as observed in PSR J0726$-$2612.

Indeed, this model is oversimplified and unlikely to apply to the real case. Whatever the mechanism responsible for the surface emission, in fact, the presence of a strong magnetic field results in some degree of anisotropy in the emitted radiation. In the case of a magnetized atmosphere, more complicated energy-dependent beaming patterns are produced: they consist of a relatively narrow pencil-beam aligned with the magnetic field, surrounded by a broader fan-beam at intermediate angles and accounting for most of the escaping radiation (see e.g. Pavlov et al. 1994). The angular pattern of the emerging intensity depends also on the local surface temperature and magnetic field, so that the morphology of the pulse profiles can be extremely variegate. Using a partially ionized hydrogen atmosphere model (Suleimanov et al. 2009) with improved opacities from Potekhin et al. (2014), we computed the expected pulse profiles, as described in Rigoselli et al. (2019). The best match with the data was obtained assuming emission from two antipodal hot spots with an effective temperature of $0.5$ MK, and $\xi=30^{\circ}$, $\chi=35^{\circ}$. In Fig. 7 we show two examples with representative values of the magnetic field, $B=4\times 10^{13}$ G and $B=6\times 10^{13}$ G. Although these pulse profiles qualitatively resemble that observed in PSR J0726$-$2612, we note that they have been computed considering only the X-ray emission from the polar caps. The addition of a contribution from an extended part of the star surface would reduce the pulsed fractions of the light curves shown in Fig. 7.

4.1 Connections with the XDINSs

Our spectral results, and in particular the presence of a broad absorption line, strengthen the similarity between PSR J0726$-$2612 and the XDINSs, for which similar spectral features have been reported (see Table 4). As it is illustrated in Fig. 8, not only the line properties, but also the best fit parameters of the continuum model are very similar to those recently reported in a systematic analysis of all the XDINS spectra with the G2BB model (Yoneyama et al. 2019).

Considerations on the age-luminosity diagram shown in Fig. 9 give even more strength to this analogy. The figure represents the bolometric luminosity of thermally emitting neutron stars against their ages, characteristic or kinematic. The luminosity of PSR J0726$-$2612 $L_{\infty}=(4.0_{-1.0}^{+4.4})\times 10^{32}$ erg s*-1* corresponds to the cold component of the G2BB fit to the phase-averaged spectrum (for $d=1$ kpc). This component is in fact representative of the cooling emission from the entire star surface (the inclusion of the hot component would not significantly change the result, adding only about 3% to the total luminosity, well within the uncertainties). The observational data for other neutron stars are displayed in Fig. 9 as in Potekhin & Chabrier (2018); most of them are taken from Viganò et al. (2013), with some updates and additions. The horizontal error bars show the uncertainties of kinematic ages, when available, otherwise the bars are replaced by arrows.

The position of PSR J0726$-$2612 in this diagram is indeed close to the group of XDINSs. Its place can be considered as intermediate between the regions occupied by ordinary neutron stars, which have either smaller luminosities or smaller ages, magnetars, which generally have larger luminosities, and XDINSs, which have somewhat smaller luminosities and larger ages. For comparison we plot two cooling curves, with heavy (nonaccreted) and light (accreted) chemical elements in the outer heat-blanketing envelope. The cooling curves are calculated for a neutron star of mass $M=1.2\,M_{\odot}$ and the dipole magnetic field inferred for PSR J0726$-$2612 ($B_{p}=6\times 10^{13}$ G) using the code of Potekhin & Chabrier (2018) with the equation of state BSk24 (Pearson et al. 2018), singlet pairing type superfluidity of neutrons and protons (according to Margueron et al. 2008 and Baldo & Schulze 2007, respectively, both in the parametrized form of Ho et al. 2015). The triplet pairing type superfluidity of neutrons is not included, because it is strongly suppressed by many-particle correlations, according to recent results of Ding et al. (2016). The latter suppression delays the onset of the Cooper pair breaking-formation mechanism of neutrino emission in the core of the neutron star and thus slows down the cooling, making the theoretical cooling curves compatible with the XDINS observations even without additional internal heating, which otherwise would be needed (e.g., Viganò et al. 2013).

While most of the XDINSs have single-peaked pulse profiles, two of them (RX J1308.6$+$2127, Hambaryan et al. 2011, and RX J0720.4$-$3125, Hambaryan et al. 2017) show double-peaked profiles similar to PSR J0726$-$2612, although with smaller pulsed fractions ($18\%$ and $11\%$, respectively). The remarkable difference between PSR J0726$-$2612 and these two XDINSs is the presence of radio emission in the former. Here we discuss the possibility that this is due an unfavourable orientation of their radio beam. Based on the radio beaming fraction of long period pulsars, Kondratiev et al. (2009) estimated that one should observe a much larger number of XDINSs ($\sim$40) to detect one with the radio beam crossing our line of sight.

We have marked in Fig. 10 the values of the angles $\xi$ and $\chi$ estimated for RX J1308.6$+$2127 and RX J0720.4$-$3125 by Hambaryan et al. (2011) and Hambaryan et al. (2017). They imply that these two pulsars are nearly orthogonal rotators ($\xi\approx 90^{\circ}$) seen with a large impact parameter $\eta=|\chi-\xi|\approx 45^{\circ}$. With the usual assumption that the radio beam coincides, or is close to, the magnetic dipole axis, such a large impact parameter can naturally account for the fact that their radio emission is not visible from the Earth. As an example, the dashed lines in Fig. 10 indicate the region where $\xi\sim\chi$ for which a radio beam with aperture of $\approx 8^{\circ}$ would be visible. Contrary to the two XDINSs, PSR J0726$-$2612 should lie inside this region. Our atmosphere model used to compute the pulse profiles of Fig 7, predicts that the radio pulse, that appears when the magnetic axis is in the plane defined by the line of sight and rotation axis, is at the phase of one of the two minima of the X-ray profile. Considering the current relative error in the radio and X-ray phase alignment (see Fig. 2), this possibility cannot be excluded.

5 Conclusions

Our analysis of XMM–Newton data of the slow, highly magnetized radio pulsar PSR J0726$-$2612 revealed the presence of a broad absorption line in its soft thermal spectrum, with parameters similar to those of the lines seen in most of the XDINSs. The X-ray pulse profile of PSR J0726$-$2612 is double-peaked and moderately energy-dependent. These findings reinforce the similarity between this radio pulsar and the XDINSs. Assuming a distance of 1 kpc, the luminosity of PSR J0726$-$2612 is $L_{\infty}=(4.0_{-1.0}^{+4.4})\times 10^{32}$ erg s*-1*. This is greater than its spin-down luminosity, as for the XDINSs (see Table 4), but it is in reasonable agreement with the expected thermal luminosity of a $\approx 200$ kyr old pulsar (see Fig. 9).

More observations are needed to reduce the uncertainty in the radio and X-ray phase alignment and better constrain the geometry of PSR J0726$-$2612. This can help to understand if the detection of radio emission in this pulsar, and not in the XDINSs with a similar double-peaked X-ray pulse profile, is due only to orientation effects.

Acknowledgements.

We are grateful to an anonymous referee for constructive suggestions. We acknowledge financial contribution from the agreement ASI-INAF n.2017-14-H.0. Part of this work has been funded using resources from the research grant “iPeska” (P.I. Andrea Possenti) funded under the INAF national call Prin-SKA/CTA approved with the Presidential Decree 70/2016. This work is based on observations obtained with XMM–Newton, an European Space Agency (ESA) science mission with instruments and contributions directly funded by ESA Member States and NASA. The work of A.Y.P. was supported by DFG and RFBR within the research project 19-52-12013. The work of V.S. was supported by Deutsche Forschungsgemeinschaft (DFG, (grant WE 1312/51-1) and by the subsidy allocated to Kazan Federal University for the state assignment in the sphere of scientific activities (3.9780.2017/8.9).