Evaluating the Evidence of Multipolar Surface Magnetic Field in PSR J0108$-$1431

TL;DR

This study investigates the surface magnetic field structure of pulsar PSR J0108-1431, finding evidence for a multipolar magnetic field based on phase offset analysis of thermal X-ray and radio emissions.

Contribution

The paper provides the first evidence for multipolar surface magnetic fields in PSR J0108-1431 through phase offset analysis of X-ray and radio emission peaks.

Findings

Bayesian analysis suggests thermal emission models are ambiguous.

Phase offset > 0.1 indicates multipolar magnetic field presence.

Thermal emission models cannot reliably estimate polar cap area.

Abstract

PSR J0108$-$1431 is an old pulsar where the X-ray emission is expected to have a thermal component from the polar cap and a non-thermal component from the magnetosphere. Although the phase-integrated spectra are fit best with a single non-thermal component modeled with a power-law (PL) of photon index $\Gamma=2.9$, the X-ray pulse profiles do show the presence of phase-separated thermal and non-thermal components. The spectrum extracted from half the rotational phase away from the X-ray peak fits well with either a single blackbody (BB) or a neutron star atmosphere (NA) model, whereas, the spectrum from the rest of the phase range is dominated by a PL. From Bayesian analysis, the estimated BB area is smaller than the expected polar cap area for a dipolar magnetic field with a probability of 86% whereas the area estimate from the NA model is larger with a probability of 80%. Due to the ambiguity in the thermal emission model, the polar cap area cannot be reliably estimated and hence cannot be used to understand the nature of the surface magnetic field. Instead, we can infer the presence of multipolar magnetic field from the misalignment between the pulsar's thermal X-ray peak and the radio emission peak. For J0108$-$1431, we estimated a phase-offset $\Delta\phi > 0.1$ between the thermal polar cap emission peak and the radio emission peak and argue that this is best explained by the presence of a multipolar surface magnetic field.

Figures

Tables

Table 1. Table 1: Pulsar J0108–1431 Parameters Summary.

Parameter	Value
Right ascension (J2000)	${01}^h{08}^{m}08\stackrel{s}{.}347016\left(88\right)$
Declination (J2000)	$-{14}^{\circ }31\prime 50\stackrel{\prime \prime }{.}1871\left(11\right)$
Position epoch (MJD)	54100
Galactic longitude/latitude ($l/b$)	$9\stackrel{\circ }{.}83$ / $-20\stackrel{\circ }{.}06$
Period ($P$)	0.807564614019(20) s
Period derivative ($\dot{P}$)	$7.704\left(12\right)\times {10}^{-17}s{s}^{-1}$
Frequency ($\nu$)	1.23829100810(3) Hz
Frequency derivative ($\dot{\nu }$)	$-1.1813\left(18\right)\times {10}^{-16}{s}^{-2}$
Epoch of timing solution (MJD)	50889
Dispersion measure (DM)	0.23 cm-3 pc
Distancea "($d$)	${210}_{-50}^{+90}$ pc
Characteristic age (${\tau }_c$)	$196$ Myr
Spin-down power ($\dot{E}$)	$4.88\times {10}^{30}$ erg $s^{-1}$
Surface magnetic field ($B_{\mathrm{surf}}$)	$2.32\times {10}^{11}$ G

Table 2. Table 2: EPIC spectral extraction parameters for PSR J0108–1431.

Detector:	EPIC-pn			MOSI	MOSII
Phase range $\varphi :$	$0-1$	$0.2-0.7$	$0.7-0.2$	$0-1$	$0-1$
Energy range (keV)	$0.15-2$	$0.15-1.2$a	$0.15-2$	$0.15-2$	$0.15-2$
Net exposure (ks)	110	55	55	126	126
Radius	$13″$	$12″$	$13″$	$12″$	$12″$
Source counts	922	309	579	218	197
Net count rate (ks-1)	$5.50\pm 0.28$	$3.87\pm 0.32$	$6.74\pm 0.44$	$1.25\pm 0.12$	$1.13\pm 0.11$
S/N ratio	20	12	15	11	10

Table 3. Table 3: Marginalized parameter distributions represented by medians with 10 − 90 10 90 10-90 percentile intervals from PSR J0108–1431 spectral fits.

	Phase-range
Parameter	$0-1$	$0.2-0.7$	$0.2-0.7$	$0.7-0.2$	$0.7-0.2$
$N_H$ (${10}^{21}$ cm-2)	${0.36}_{-0.12}^{+0.15}$	${1.46}_{-0.64}^{+1.15}$	${0.21}_{-0.12}^{+0.17}$		${0.022}_{-0.019}^{+0.052}$
$\Gamma$	${2.93}_{-0.23}^{+0.25}$	${5.0}_{-1.0}^{+1.4}$	—	${2.45}_{-0.27}^{+0.33}$	—
PL norm ($N_{-6}$)a	${1.62}_{-0.42}^{+0.51}$	${1.51}_{-0.60}^{+0.83}$	—	${3.55}_{-0.42}^{+0.44}$	—
$F_{0.3-7keV}^{\mathrm{unabs}.}$b	${1.41}_{-0.12}^{+0.15}$	${3.1}_{-1.6}^{+8.1}$	—	${1.67}_{-0.20}^{+0.24}$	—
$k{T}_{BB}$ (eV)	—	—	${115}_{-12}^{+12}$	—	${164}_{-13}^{+14}$
BB normc	—	—	${5.6}_{-2.3}^{+4.9}$	—	${1.48}_{-0.41}^{+0.60}$
BB area ratiod	—	—	${0.38}_{-0.16}^{+0.34}$	—	${0.10}_{-0.03}^{+0.04}$
$L_{BB,bol}$ (${10}^{28}$ erg s-1)	—	—	${5.62}_{-0.98}^{+1.59}$	—	${6.16}_{-0.52}^{+0.57}$
${\chi }_{\nu }^2$	1.2	1.5	1.2e		1.7

Table 4. Table 4: Timing extraction parameters.

Energy Range	GTI	Radius	(Net/Total)	S/N	$Z_2^2$
(eV)	(ks)
$150-300$	127	$22\stackrel{\prime \prime }{.}5$	244 / 413	11.6
$300-400$	127	$22\stackrel{\prime \prime }{.}5$	122 / 190	8.6
$400-500$	127	$10\stackrel{\prime \prime }{.}5$	81 / 93	8.3
$500-600$	127	$10\stackrel{\prime \prime }{.}5$	61 / 72	7.2
$600-700$	118	$15\stackrel{\prime \prime }{.}5$	54 / 67	6.6
$700-900$	127	$12\stackrel{\prime \prime }{.}5$	67 / 95	6.8
$900-2000$	127	$7\stackrel{\prime \prime }{.}5$	76 / 126	6.8
$150-300$	92	$7\stackrel{\prime \prime }{.}5$	83 / 95	8.5	28.8
$300-500$	127	$10″$	146 / 172	11.1	17.6
$500-600$	84	$5″$	22 / 23	4.6	24.9
$600-1000$	109	$22\stackrel{\prime \prime }{.}5$	131 / 226	8.5	12.8
$1000-2000$	127	$20″$	51 / 144	4.1	36.2

Table 5. Table 5: Timing Solution for PSR J0108–1431.

Method	EPOCH	$F0$	$F1$	$F2$	JUMP
	(MJD)	(Hz)	($\times {10}^{-16}$ s-2)	($\times {10}^{-27}$ s-3)	(s)
Original	50889	1.23829100810(3)	-1.1813(18)	-8.9(77)	—
Expecteda	55728	1.23829095793322	-1.2185	-8.9	—
Fit Radiob	55728	1.2382909586474(77)	-1.1789(48)	-43(26)	—
Fit X-rayc	55728	1.238290958651(12)	-1.1822(67)	-53(38)	-0.03995(29)

Equations

$$A_{\mathrm{d},\mathrm{pc}}=2{\pi }^2{R}_{\mathrm{NS}}^3/(cP)$$

$$\begin{array}{rl}P(d)=& H(D-d)\exp \left[-0.5\times {\left(\frac{d-D}{D_{\min }}\right)}^2\right]\\ +& H(d-D)\exp \left[-0.5\times {\left(\frac{d-D}{D_{\max }}\right)}^2\right],\end{array}$$

$$H(x)=$$

$$\sin (\alpha )/\sin (\beta )=\mid d\Psi /d\Phi {\mid }_{max}$$

$${\sin }^2(\rho /2)=\sin (\alpha +\beta )\sin (\alpha ){\sin }^2(\varphi /4)+{\sin }^2(\beta /2)$$

$$f(x)=A0+A\times \sin [2\pi (x-{\varphi }_0)]$$

Full text

††pubyear: 2019††pagerange: Evaluating the Evidence of Multipolar Surface Magnetic Field in PSR J0108–1431–Evaluating the Evidence of Multipolar Surface Magnetic Field in PSR J0108–1431

1 Introduction

Rotation-powered pulsars (RPPs) are rapidly rotating, highly magnetized neutron stars (NSs) that have active magnetospheres, where efficient plasma production and acceleration processes can occur. Due to the rapid rotation, the magnetosphere is divided into regions of open and closed magnetic field lines depending on whether or not they extend beyond the light cylinder — an imaginary cylinder around the rotation axis where the corotation speed equals the speed of light. In the standard model of pulsars, a strong co-rotating electric field populates the NS exterior by pulling out charges from the NS surface (Goldreich & Julian, 1969) and/or creating electron-positron pairs through magnetic pair production (Sturrock, 1971). The charged particles achieve a force-free (FF) configuration due to the strong electromagnetic forces that vastly exceed the inertial and dissipative forces. Under such conditions, the minimum charge needed to screen the electric field is the Goldreich-Julian number density given by $n_{\rm GJ}=\Omega\cdot B/2\,\pi\,e\,c$, where, $B$ is the magnetic field, $\Omega=2\,\pi/P$ ($P$ is pulsar period), $e$ is the electric charge, and $c$ is the speed of light in vacuum. The closed field line region is filled with dense plasma which screens the electric field and co-rotates with the star. The FF conditions also provide prescriptions for the charge density and current density in the open field line regions. However, the FF conditions need to be relaxed in the open field line regions to explain observations of energetic particle winds and high-energy radiation that are powered by the pulsar’s spin-down.

The evidence for particle wind comes from observations of bright wind nebulae around young pulsars (Kargaltsev et al., 2015), and modeling such wind nebulae requires high pair multiplicity factors, $\kappa\lesssim 10^{5}n_{\rm GJ}$ as suggested by the pair cascade models (see e.g. de Jager 2007; Timokhin & Harding 2019). The magnetospheric origin of the electromagnetic radiation is apparent from the pulsed and co-rotating radiation patterns observed from the radio to the $\gamma$-rays in several pulsars. The open field lines regions, hence, require “gap” regions where copious pair creation, acceleration, and cascade takes place, which leads to a highly-relativistic flow of charges as particle winds and the electromagnetic broadband radiation. Although the modeling of FF magnetosphere has been successful in providing global solutions for the magnetic field structure and the plasma dynamics in the pulsar magnetosphere (see for e.g. Michel & Li 1999; Contopoulos et al. 1999; Spitkovsky 2011; Pétri 2016), it is yet to incorporate self-consistent particle acceleration and radiation models (Cerutti & Beloborodov 2017; and references therein). Observations that constrain the pair cascade gap regions and the local magnetic field geometry are crucial to determine regions where the FF conditions are broken.

To obtain such observational constraints on the magnetic field structure, in this work we will focus on emission from normal RPPs with periods longer than 100 ms. Observations of a few pulsars in this category show both incoherent X-ray and coherent radio emission. The X-ray spectra are usually modeled to have a combination of thermal and non-thermal emission components. The thermal emission is thought to originate from the polar cap111Younger pulsars, with characteristic age $\tau\lesssim 1$ Myr, have an additional X-ray emission component from the entire stellar surface as they lose their heat of formation.. The physical location of the non-thermal emission, on the other hand, is uncertain and various models suggest that the emission can originate from different parts of the magnetosphere (Becker 2009; Cerutti & Beloborodov 2017). The coherent radio emission has been constrained by observations to originate from regions below 10% of the light cylinder (see for e.g. Mitra (2017) and references therein). These features of RPP emission are best explained by Ruderman & Sutherland (1975), RS model of pulsar emission, where it was postulated that a charge-starved “inner vacuum gap” (IVG) region exists just above the pulsar polar cap (PC) where pair plasma can be created, accelerated, and cascaded. The efficiency of pair creation depends on the curvature of the magnetic field lines in the IVG, and to get $\kappa\sim 10^{5}$, the dipolar magnetic field radius of curvature is too large and hence inadequate (Timokhin & Harding, 2019). Instead, the presence of strong multipolar magnetic field lines with smaller radii of curvature is more efficient at increasing the pair multiplicity.

The electric field in the IVG accelerates the charged particles of one sign upwards along the open-field lines and the oppositely charged ones downwards onto the polar cap (Ruderman & Sutherland, 1975). The particle bombardment on the polar cap heats up the surface to temperatures $kT\gtrsim 0.1$ keV (Ögelman 1991; Cheng & Ruderman 1980; Harding & Muslimov 2001). This has been observed as thermal emission from emitting regions of sizes much smaller than the NS radius (see e.g. Pavlov et al. 2002; Geppert 2017; Rigoselli & Mereghetti 2018 and references therein). This emission flux depends on the energy dumped by the accelerated charged particles and the PC area. For a global dipolar magnetic field, the PC area is calculated by finding the locus of the last closed field lines on the NS surface. A multipolar magnetic field geometry has dipolar and other higher order contributions near the surface while maintaining a purely dipolar structure away from the surface (Gil et al., 2002). Due to magnetic flux conservation, such multipolar magnetic field structures are expected to have smaller PC area than the PC area for a pure dipole. Since we detect X-ray thermal emission from localized hot-spots (presumably the PCs) in some RPPs, constraining their emission size can provide evidence for/against a purely dipolar field.

In the case of a blackbody, the emission is isotropic and the flux is directly proportional to the emitting area. The isotropic nature of the BB at all energies produces pulse profiles independent of energy and dependent only on the geometry and inclination of the emitting region. Hence, measures of profile shape, such as pulsed fraction – the ratio of counts in the pulse over the total counts, are expected to remain constant over energy for a BB. Once the BB nature of the PC emission is established, its area is estimated from fitting a BB model component to the observed spectrum or by modeling the pulse profile using geometric modulation of a hot spot on a rotating sphere (Pechenick et al. 1983; Page 1995). The emission area can then, in principle, be compared with the theoretical PC area for a magnetic dipole to establish the presence of a multipolar magnetic field geometry.

In reality, however, the emission from the PC can deviate from a BB due to a variety of effects. The NS crust could have anisotropy in thermal conductivity due to the strong magnetic fields that constrain the motion of the electrons (Pons et al. 2009; Reisenegger 2009, Geppert et al. 2004). The star could have a thin atmosphere on the surface, containing elements with energy levels modified by the strong magnetic field, which reprocesses the BB emission (Zavlin & Pavlov 2002; Pavlov & Shibanov 1978). The strong backflow of particles can modify the thermal emission through inverse Compton scattering (Kardashëv et al. 1984, Sturner et al. 1995, Baring et al. 2011). In such cases, determining the emission process and estimating the polar cap area are not straightforward and the solutions may not be unique. In this paper we will often single out BB from thermal emission models where the former is an un-modified blackbody emission whereas the latter, thermal emission, can be either a BB or BB modified by one of the above-mentioned processes.

The PC area estimation procedure outlined above to check for multipolar magnetic fields requires definite proof that the PC emission is a BB. However, there is no definitive proof uniquely identifying the thermal emission from pulsars as BB. Phase-integrated spectra show thermal emission as one of its components which are routinely fit with a BB model. But, with current high signal-to-noise (S/N) spectra, the fits from neutron star atmosphere models are statistically indistinguishable from a BB fit (Bogdanov, 2013). Rather, pulse profiles in the energy range where thermal emission is dominant are often show large pulsed fractions and modeling shows deviations from the BB profiles (Bogdanov, 2013). Due to low statistics of most X-ray detected RPPs, tests of isotropic emission (non-evolution of the pulsed fraction with energy) are inconclusive (e.g., Hermsen et al. 2018). Hence, we have not yet identified a single RPP with definite proof of BB emission component from its PC.

Our goal is to obtain robust evidence for multipolar magnetic fields by performing a uniform analysis of thermal emission from non-recycled RPPs. There have been attempts to find evidence of multipolar magnetic fields by estimating the PC areas (Geppert 2017; Rigoselli & Mereghetti 2018). However, these works rely on a diverse body of literature, which suffers from one or more of a number of problems (See Section 7 for a detailed critique). Foremost of all, convincing evidence for BB emission, even in some cases the presence of thermal components, is not provided. Further, the confidence limits on the BB area obtained from spectral fitting are often underestimated due to the fixing of secondary model components or some of the parameters such as distance and NS radius. We would like to address these issues by attempting a uniform statistical analysis of all thermally emitting, non-recycled RPPs, and determine whether unambiguous proof of BB PC emission can be obtained with the currently available data. The analyses and tests will be customized to characterize the thermal emission from RPPs and estimate the emitting region areas where possible.

In cases where the thermal emission cannot be uniquely modeled, we will find evidence for/against a purely dipolar magnetic field by estimating the offsets between the thermal X-ray and radio peaks (X-R offsets). The constraint on the radio emission heights ($\lesssim 10\%$ of light-cylinder) and the success of the rotating vector model (RVM; Radhakrishnan & Cooke 1969), for slowly rotating RPPs, is consistent with a dipolar magnetic field structure at the radio emission region. A star-centered, global dipolar magnetic field geometry will hence show X-ray thermal PC emission peak coincident with the radio peak, after accounting for the effects of aberration and retardation. Alternatively, measurements of significant X-R offsets will provide qualitative evidence for multipolar magnetic field structure close to the NS surface (see e.g., Gil et al. 2002; Melikidze et al. 2009; Szary et al. 2015).

To begin with, in this work, we study the pulsar PSR J0108–1431 (J0108 hereafter), an old RPP where an earlier analysis had shown the possibility of a thermal emission component in addition to a dominant non-thermal component (Posselt et al. 2012; BP2012 hereafter). It is representative of a class of old and cooled NSs with characteristic ages $\tau_{c}\gtrsim 1$ Myr (e.g., PSR B1929+10; Becker et al. 2006, PSR B0950+08; Pavlov et al. 2017, and PSR B1133+16; Szary et al. 2017), where the thermal emission from the bulk of the surface is not seen in the X-ray band. The overall X-ray emission (thermal PC and non-thermal) gets fainter with age as the spin-down power reduces (Li et al., 2008). Much of the analyses in these systems are restricted to obtaining the best S/N pulsations and spectra and testing simple emission models with few parameters. Through spectral and timing analysis, we attempt to find stronger evidence for thermal PC emission from the pulsar and obtain conservative estimates for the model parameters. Then we use the thermal PC emission area and phase of peak flux to evaluate the evidence for a multipolar surface magnetic field.

PSR J0108–1431 (Tauris et al., 1994) is a slow-rotating ($P=0.808$ s) pulsar with a characteristic age $\tau_{c}=196$ Myr, at a parallax measured distance of $210_{-50}^{+90}$ pc (Table 1; Deller et al. 2009; Verbiest et al. 2012). It has a spin-down power ($\dot{E}$) of $5\times 10^{30}$ erg s*-1* and an estimated surface magnetic field ($B_{\rm surf}$) of $2.3\times 10^{11}$ erg s*-1*, assuming the power loss is entirely via magnetic dipole radiation. It was detected in the X-rays with Chandra (Pavlov et al., 2009) and later re-observed with XMM-Newton (BP2012). The pulsar showed a relatively soft phase-integrated spectrum well-fit with a PL of photon index $\Gamma\approx 3$ in both these data. BP2012 showed that a BB component could be fit only if the PL photon index is fixed at 2. Although X-ray pulsations were detected in the XMM-Newton data, a phase-resolved spectral analysis was not attempted due to low counts.

We re-analyze the data using improved energy-dependent events extraction to obtain pulse profiles with higher pulsations significance. This has allowed us to select the optimal energy range where the thermal emission can be better distinguished from the non-thermal emission. Hence, we were able to detect the thermal emission with higher significance and determine reliably the phase at which the emission peaks.

2 Observations and Data Analysis

We use the archival XMM-Newton data of PSR J0108–1431 from June 15, 2011 (MJD 55727). The European Photon Imaging Cameras (EPIC) were operated in full-frame mode during the continuous 127 ks exposure. We mainly work with the EPIC-pn (Struder et al. 2001) data which has a 73.4 ms (less than pulsar period $P=0.808$ s) time-resolution for timing and phase-resolved analyses. The EPIC-MOS detectors, operated in full-frame mode, have a time-resolution of 2.6 s (greater than the pulsar period) and hence, were used only for phase-integrated spectral analysis. The data reduction was performed with XMM-Newton Science Analysis Software (SAS) ver. 16.1.0 (Gabriel et al. 2004), applying standard data reduction and events filtering and extraction tasks.

We identify time-intervals that are affected by soft proton flaring background by extracting a count rate light curve of single pixel events (from the entire CCD) with energies greater than 10 keV. At these high energies, the individual source contributions are significantly reduced or negligible and the flaring background contribution dominates. Soft proton flaring background had affected $\sim 30\%$ of the exposure with count rates exceeding the quiescent level by over $3\sigma$ and approximately $5\%$ with count rates exceeding the quiescent level by $10\sigma$ (Figure 1). The good time intervals (GTIs) are obtained by filtering out time ranges when the flaring count rate exceeds a chosen threshold.

The source is detected only in the $0.15-2$ keV energy range with weak detection between $1-2$ keV and no significant contribution above 2 keV (Figure 2). The sensitivity of the EPIC-pn detector is energy dependent and for a typical pulsar spectrum, the count rates decrease steeply with energy. This leads to non-uniform signal-to-noise (S/N) with energy for events extracted from a constant size region. In order to obtain a high S/N events extraction for timing and spectral analysis, we optimize events extraction over the three parameters – energy range, GTIs, and extraction radius.

To extract an optimal number of photon events at all energies, we determine the S/N independently for the energy bins: $150-300$ eV (first bin), 100 eV bins in the $300-1000$ eV range, 1 keV bins in the $1-7$ keV range, and a $7-10$ keV bin (last bin). This selection of energy bins is arbitrary, but follows the condition that the bin widths exceed the detector’s minimum energy resolution (by a factor $\geq 2$ below 1 keV and $\sim 7-20$ above 1 keV). Large bins at high energies are preferred due to rapidly decreasing source counts with energy. In these energy ranges, we determine S/N values over a range of GTI selections and extraction region radii. The $55\arcsec$ radius circular region used for background events extraction is shown in Figure 2. Ten sets of GTIs are obtained by choosing flaring count-rate cut-offs between the $3\sigma$ limit over the quiescent level and the highest count rate (Figure 1) while ensuring a linear increase in total exposure time. The radii of source extraction apertures range from $5\arcsec$ to $22\aas@@fstack{\prime\prime}5$ in increments of $2\aas@@fstack{\prime\prime}5$. The center of the extraction circle is chosen by taking the median physical X and Y positions from the events extracted from an initially provided $22\aas@@fstack{\prime\prime}5$ radius extraction region centered approximately on the source.

We obtain optimal parameters for the highest S/N source events extraction from the matrix of S/N values calculated over the range of GTIs and extraction radii for each energy bin. For J0108, the method did not benefit the spectral analyses and hence the spectra were extracted following the standard procedure using a single GTI selection and extraction region over the full energy range. But, as shown below, we achieved high significance pulsations not only in the full $0.15-2$ keV range but also in the two energy ranges below and above 0.7 keV. As a result, we obtained reliable profiles and this prompted further phase-resolved spectral analyses.

3 Phase-Integrated Emission

X-ray emission from neutron stars can originate from a few different sites:

(i) The heat of formation and any subsequent internal re-heating is lost via thermal radiation from the bulk of the stellar surface (Figure 3a). NS thermal evolution models predict cooling times of the order of 1 Myr beyond which the surface temperature is too low to produce significant X-rays (Gnedin et al. 2001; Yakovlev et al. 2005; Page et al. 2006). Thermal soft X-rays from regions of sizes comparable to the NS radius have been observed from young pulsars (e.g., Geminga; Halpern & Ruderman 1993, PSR B1055–52; Oegelman & Finley 1993, PSR B0656+14; Possenti et al. 1996, Vela; Pavlov et al. 2001).

(ii) The polar cap region on the NS surface is heated by highly-relativistic charged particles that are accelerated towards the surface (Figures 3a,3b). X-ray detected pulsars, both young and old, often show a relatively hotter and localized thermal component in their spectra (e.g., see De Luca et al. 2005). This apparent hot spot is often identified with the polar cap. The polar cap heating continues as long as the NS can accelerate particles in the inner vacuum gap (Ögelman 1991; Cheng & Ruderman 1980).

(iii) Charged particles accelerated in regions/gaps inside the open field lines can emit X-rays via synchrotron Zhang et al. 1998 or possibly inverse-Compton radiation (Lyutikov, 2013). Alternatively, certain regions within the return current can also accelerate charged particles and produce X-ray emission (Cerutti & Beloborodov 2017; and references therein). Non-thermal emission is very commonly observed in the X-ray spectra of RPP (Li et al. 2008, and references therein).

A phase-integrated spectrum is obtained by using source events from all detected rotational phases. The individual emission components are better resolved when each of them dominate at a different energy range. Figure 3a-i illustrates a case wherein, for a young pulsar, the spectrum is dominated by the thermal emission from the bulk surface at low energies ($\lesssim 0.5$ keV), the PC hot spot emission in the intermediate energies ($\approx 0.5-1.5$ keV), and the non-thermal emission at higher energies ($\gtrsim 1.5$ keV). For an old pulsar, an ideal scenario (Figure 3b-i) will have the thermal hot-spot emission dominating the lower energies ($\lesssim 1$ keV) and the non-thermal component dominating the higher energies ($\gtrsim 1$ keV). For the purposes of this illustration, we use BB to represent the thermal component and PL for the non-thermal.

The spin-down power of pulsars reduces by as much as $\sim 3-4$ orders of magnitude over the $\sim 10^{7}$ years decrease in characteristic age, as seen from the $P-\dot{P}$ diagram of pulsars. Additionally, as the pulsar slows down, the light-cylinder radius increases resulting in reduced polar cap area. As a result, the emission from old pulsars is considerably weaker than those of younger pulsars for comparable distances, assuming no increase in PC heating. As explained earlier, observed thermal emission area can be considerably smaller than the PC area for a dipole due to the orientation/geometry of the pulsar and the presence of multipolar magnetic fields. This could lead to scenarios where the thermal component in the spectra dominates only in a very narrow energy range (Figure 3b-ii) or is simply too weak compared to the non-thermal component (Figure 3b-iii). In such cases, with low S/N data as is typically the case, the phase-integrated spectrum cannot differentiate between a two-component thermal+non-thermal emission model and a single non-thermal emission model (e.g., J0108 and the set of old NSs listed above).

The source, J0108, is detected only in the $0.15-2$ keV range, suggesting a soft spectrum. The EPIC-pn and MOS spectral extraction parameters optimized over GTIs and extraction aperture are shown in Table 2. We group the spectra with at least 1 count per energy bin while ensuring that the bin width is at least 1/3 of the detector resolution full width at half maximum (FWHM).

We perform spectral fitting using XSPEC 12.10.0c. We model the spectrum with a simple power-law (PL; powerlaw in XSPEC), modeling the interstellar absorption using the Tübingen-Boulder model (Tbabs in XSPEC; Wilms et al. 2000) with solar abundance table wilm from Wilms et al. (2000), and photoelectric cross-section table bcmc (Balucinska-Church & McCammon 1992; Yan et al. 1998).

We perform Bayesian fits to the spectra using the statistic lstat222https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html in XSPEC which constructs the likelihood function by combining the source and background counts (without subtraction) assuming Poisson distribution for both. The posterior probability distribution is marginalized over the background rate parameter to retain only the distributions of the emission model parameters. We select the Goodman-Weare algorithm (Goodman & Weare, 2010) and set up Monte Carlo Markov chain (MCMC) with 50 walkers taking 2000 steps each (1000 steps for burn-in) to sample the posterior probability distribution over the model parameter space. The PL fit is similar to the one obtained in BP2012. A $\Gamma=3$ PL fits the spectra well with no significant systematics in the residuals. In Figure 4, we show the spectral fit with residuals and in Figure 5, we show the marginalized individual and joint posterior probability distributions of the fit parameters.

Other single-component models such as blackbody (BB: bbodyrad in XSPEC) and neutron star atmosphere (NSA: nsa or NSMAXG: nsmaxg models in XSPEC) fit worse than PL, with high residuals in the high energy bins ($\gtrsim 1$ keV). Additional components to the PL are not required to explain the phase-integrated spectrum, due to the low S/N of the data. However, an additional thermal component (BB) is allowed with a very low probability and seen to dominate the emission at energies $\lesssim 0.8$ keV. We conclude that the phase-integrated spectra of J0108, due to their low S/N and possibly a relatively weaker thermal component, cannot confirm or reject the canonical model of two-component emission from old pulsars.

4 X-ray pulsations from the pulsar

The 73.4 ms time-resolution of the EPIC-pn detector allows us to assign rotational phase of photon events with an uncertainty $\Delta\phi\approx 0.09$. Following the events extraction procedure explained in section 2, we extracted 1057 counts from the source region in the $0.15-2$ keV energy range with different GTIs and extraction radii for each energy bin (Table 4). This ensures high S/N events extraction in the selected energy bins. The events were barycentered using the SAS task barycen. We use the $Z_{n}^{2}$ test (for up to n=4 harmonics) to find the pulsation frequency within a 0.002 Hz window around the pulsar frequency $\nu_{\rm exp}=1.2382909579(8)$ Hz, extrapolated from the known ephemeris at MJD 50889, where, the $1\sigma$ uncertainty in the last digit is given in parenthesis. We detect pulsations at $\nu=1.23829043(53)$ with $Z_{2}^{2}=30.2$ ($4.4\sigma$) and statistically significant contributions up to 2 harmonics. The detected frequency is within $1.5\sigma$ of the expected frequency. For the observation span of 127 ks, $Z_{n}^{2}$ is only sensitive to frequency derivative $\gtrsim 8\times 10^{-6}$ Hz s*-1* and frequency second derivative $\gtrsim 8\times 10^{-16}$ Hz s*-2*, both of which are orders of magnitude greater than the pulsar’s frequency derivative and second derivative (see Table 1).

In order to increase the pulsations significance, we repeated the source events extraction over the GTIs, energy ranges, and extraction aperture sizes by maximizing the $Z_{2}^{2}$ statistics instead of the conventional S/N maximization. We extracted 660 counts from the source region in the $0.15-2$ keV range with varying GTIs for different energy bins, and extraction aperture radii varying between $5\arcsec-15\arcsec$. These events show significantly stronger pulsations at $\nu=1.23829044(25)$ with $Z_{4}^{2}=85.2$ ($7.8\sigma$) and statistically significant contributions up to 4 harmonics (Figure 6). The $Z_{2}^{2}=71.7$ is also significantly higher than the value obtained from simple S/N optimization.

Phase-folding the source and background events with the frequency of the detected pulsations, we obtain a binned pulse profile choosing MJD 55727.93394840 (XMM reference MJD 50814.0 + mission elapsed time 424563893.142/86400 days) as phase zero (Figure 6). The frame-time of the detector mode spans $\approx 0.09$ in phase and hence, for a profile with 10 bins, the uncertainty in phase for any event is close to the bin width. To obtain a reliable binning-independent pulse profile, we use Kernel Density Estimation (KDE; Feigelson & Babu 2012) smoothing using rectangular (boxcar function) kernels for each event with the edges matching the event frame-times. This assigns a uniform probability of an event detection within its frame-time. Such probability distributions for all events are then superimposed to obtain the probability density function (PDF) of event detection at any phase. This PDF is multiplied with the average counts per bin of any binning scheme to obtain the corresponding smoothed profile.

The binned and smoothed phase-integrated profile in the 0.15-2 keV range is shown in Figure 6. The profile is strongly peaked around phase $\phi=0$ for our chosen folding time reference. Such narrow-peaked (non-sinusoidal) profiles are characteristic of non-thermal emission or highly beamed thermal emission. BB emission, on the other hand, shows broad sine-like pulse with large duty cycle when special and general relativistic effects are included. We attempted to detect thermal pulsations assuming the presence of a thermal component in the phase-integrated spectra. In trial spectral fitting with combined PL+BB model, the thermal component dominates at energies $\lesssim 0.8$ keV. Hence, we produced energy-resolved pulse profiles in the $0.15-0.7$ keV, and $0.7-2$ keV ranges (Figure 7).

We estimate the phase of pulse profiles and its uncertainties using Bootstrap resampling (Feigelson & Babu, 2012), assuming Poisson distributed counts in each bin and taking the observed values to be the mean. The peak phase estimate is shown as Box plots in Figures 6b and 7), where the bars on the box represent the three quartiles (25%, 50%, and 74%), the whiskers connected to the box at $1.5\times$ the inter-quartile range represent the probable limits of the distribution and the points beyond are the outliers. The overall pulse peak is maintained at $\phi\approx 0$ in the selected energy ranges with no statistically significant shift in peak phase.

In addition to the main pulse seen at all energy ranges near phase 0, the profiles show a possible second pulse around phase 0.5 which is dominant only in the $0.15-0.7$ keV range. This hints at the possibility of a separate soft-emission component. We assess the significance of the overall profile change in the two energy ranges by applying the Anderson-Darling test (Anderson & Darling, 1954) to the binned profiles. The null hypothesis that profiles are samples obtained from the same underlying distribution is not rejected, possibly due to the low number of counts. Hence, the change in profile from double-peaked in $0.15-0.7$ keV to single-peaked between $0.7-2$ keV is not sufficient to infer two phase-separated emission components. However, we use spectra extracted from the two phase ranges, $0.2-0.7$ and $0.7-0.2$ to further explore the emission components.

5 Phase-separated Spectral Components

From the pulse profile shapes, we inferred soft emission in the $0.2-0.7$ phase range and relatively harder emission in the rest of the phases ($0.7-0.2$). We extract spectra in these two broad phase ranges from the EPIC-pn detector. A finer phase-resolved spectral extraction is not reliable due to the low number of source counts and the coarse time-resolution. The spectrum extracted from the $0.2-0.7$ phase range (S1) has $213\pm 18$ net source counts in the $0.15-1.2$ keV range and none over 1.2 keV (Table 2). The spectrum extracted from the $0.7-0.2$ phase range (S2), on the other hand, has $371\pm 24$ net source counts in the $0.15-2$ keV range.

5.1 Blackbody Emission Model

Single-component models, either thermal or non-thermal, are sufficient to explain the spectra. We obtained the best fit (Figure 8) and constraints (Table 3; Figure 9) on the model parameters when simultaneously fitting a BB for S1 and a PL for S2, while using a shared ISM absorption parameter ($N_{\rm H}$ tied). The resulting constraint on $N_{\rm H}$ is stronger due to the improved statistics from using data from both the spectra than otherwise obtained from individual spectral fits. This also better constrains the emission model parameters. A PL fit to S1 produces unrealistically steep photon index $\Gamma\approx 5$ and results in a bad overall fit ($\chi_{\nu}^{2}=1.5$), while a BB fit to S2 is overall bad ($\chi_{\nu}^{2}=1.7$) with large systematic residuals.

BB emission flux is directly proportional to the emitting area and the spectral shape is uniquely determined by the temperature parameter. The polar cap area on a NS of radius $R_{\rm NS}$ and spin period $P$ for an spin-aligned magnetic dipole is

[TABLE]

where, $c$ is the speed of light (Sturrock, 1971). For J0108, assuming $R_{\rm NS}=10$ km and mass of the NS $M_{\rm NS}=1.4\,M_{\sun}$, $A_{\rm d,pc}=0.0815\,{\rm km}^{2}$, which for an observer at infinity is $A_{\rm d,pc}^{\infty}=A_{\rm d,pc}/g_{\rm r}^{2}=0.139\,{\rm km}^{2}$, where, $g_{\rm r}=\sqrt{1-2.952\,M_{\rm NS}/R_{\rm NS}}$ is the gravitational redshift parameter. Through Bayesian analysis, we obtain a representative sample of the probability distribution function (PDF) for the emitting area (the marginalized posterior probability distribution conditional on the observed spectrum) assuming the nominal distance estimate of 210 pc. The median value $A_{\rm d,pc}^{\infty}=0.032_{-0.014}^{+0.031}$ km2, with uncertainties obtained from the 10-90 percentile interval. We estimate the probability that the measured BB area is lower than the nominal value of $A_{\rm d,pc}^{\infty}$, $P(A<A_{\rm d,pc}^{\infty})=99.8$%. We propagate the pulsar’s rather large and asymmetric distance uncertainties, used in the BB area estimates, by sampling a distribution of the form

[TABLE]

where, $D_{\rm min},D,D_{\rm max}$ are 160, 210, and 310 pc, respectively (Verbiest et al., 2012). The probability that the measured BB area is lower than $A_{\rm d,pc}^{\infty}$, after including the distance uncertainty, is $P(A<A_{\rm d,pc}^{\infty})=86$%. The ratio of the measured hot spot area to $A_{\rm d,pc}^{\infty}$ is $0.48_{-0.30}^{+0.78}$, with uncertainties quoted at 90% probable limits about the median value of the distribution.

We tried to constrain the phase of the thermal peak by comparing thermal emission component contribution in the two spectra S1 and S2. We simultaneously fit a BB model to S1 and a PL+BB model to S2 while keeping the parameters $n_{H}$ and $kT$ tied-up between the two spectra (Table 3). The area parameter estimates in the two phase ranges show larger area values for spectrum S1 compared to S2 (Figure 8). We calculate a 98% probability that the BB area is larger in the $0.2-0.7$ phase range when compared to that in the $0.7-0.2$ phase range.

5.2 Beamed Thermal Emission Models

A BB is not the only thermal emission model that can explain the neutron star surface radiation. Local magnetic field structures can produce beamed radiation by introducing anisotropy in the conductivity of the stellar crust, radiation transfer in a layer of NS atmosphere, or charged particle motion that Comptonize the radiation. To compare the effect of beaming on observed emission areas, we fit neutron star atmosphere models (NSA and NSMAXG) to the soft spectrum extracted in the $0.2-0.7$ phase range.

We fit for the effective temperature and the normalization parameters in the NSA model while freezing the NS mass at $1.4\,M_{\sun}$, NS radius at 10 km, and the surface magnetic field strength at $10^{12}$ G. The overall fit of the PL+NSA model to the spectra (NSA to S1 and PL to S2) is comparable to the PL+BB fit ($\chi_{\nu}^{2}=1.2$). The $10-90$% parameter range about the median for the effective temperature $T_{\rm eff}=9.0_{-1.6}^{+2.1}\times 10^{5}$ K and the emission area to NS area ratio $A_{\rm ratio}=2.6_{-1.7}^{+5.3}\times 10^{-4}$ ($A_{\rm pc}/A_{\rm NS}=6.3\times 10^{-5}$). Comparing the fit parameters — temperature and emission area — between the two models, we see that the spectrum is unable to differentiate between a set of BB fits and an equivalent set of NSA fits with colder temperatures and larger emitting areas (Figure 10). The comparison shows that spectra from a beamed, reprocessed thermal emission model such as from a NS atmosphere, cannot be statistically distinguished from a BB. In this instance, fitting BB to a thermal spectrum processed by a NS atmosphere would lead to under-estimated emission areas and over-estimated temperatures.

5.3 Absorption-like features in the thermal component

The spectrum in the phase-range $0.2-0.7$ shows an absorption-like feature at approximately 0.3 keV. The feature, seen as un-modeled residuals in thermal BB or NSA model fits, is hard to explain with purely instrumental systematics. The $\sim 40\pm 20$% residuals are significantly greater than the $3\sigma$ upper-limit on the systematic uncertainty in the effective area of $4\%$ estimated for the EPIC-pn detector333http://xmm2.esac.esa.int/docs/documents/CAL-TN-0018.pdf. Moreover, the feature is only seen in the spectrum from one-half of the rotational phase (Figure 8).

We refer only to the model components fit to spectrum S1 below, even though we always simultaneously fit a PL to S2. Introducing an absorption line modeled with a Gaussian profile (Gabs in XSPEC) to the NSA model, we see a marginal improvement in the test statistic ($\chi_{\nu}^{2}=1.1$ for Gabs$\times$NSA versus $\chi_{\nu}^{2}=1.2$ for NSA). We estimate the central energy of absorption $E_{\rm c}=0.33_{-0.08}^{+0.05}$ keV with equivalent width $EW=0.105\pm 0.05$ for NSA continuum with $\log T_{\rm eff}=5.92_{-0.07}^{+0.05}$ and area ratio $A_{\rm ratio}=4.1_{-1.8}^{+6.8}\times 10^{-4}$.

Note on statistical model selection:

The models (BB/NSA/NSMAXG)+PL and Gabs(NSA+PL), fit simultaneously to the phase-separated spectra, all produce statistically acceptable fits. A measure such as the Deviance Information Criteria (DIC; Spiegelhalter et al. 2002; Gelman et al. 2003) fails to statistically distinguish between the models, with the values lying in a narrow range of 64–71444A difference of $\geq 10$ is typically considered strong evidence and $5-10$ may be considered as substantial evidence in favor of the model with the lower DIC..

6 Thermal X-ray and Radio Phase Alignments

In the earlier sections, we have demonstrated that thermal emission mechanism cannot be uniquely identified as BB and consequently, the PC area can be very different depending on the emission model assumed. Hence, X-ray spectra cannot unambiguously provide evidence for the surface multipolar magnetic field in J0108.

An alternative way of accessing the presence of multipolar magnetic fields is to look for phase offsets between the time-aligned radio and X-ray pulse profile. The radio emission is known to arise a few hundred km above the neutron star surface. Hence, if the thermal X-ray emission originates from regions on the surface with a multipolar magnetic field, the higher order multipoles fall off rapidly with height, and at the radio emission height, the magnetic field is purely dipolar. Thus, in the presence of a multipolar magnetic field, the time-aligned thermal X-ray and radio profiles can show significant phase offsets. In this section, for J0108, we first check the validity of the assumption that the radio emission arise from regions where the magnetic field is dipolar and then time-align the radio profile with the thermal X-ray peak to obtain the phase-offset.

6.1 Radio Emission, profile, geometry, and Emission heights

Radio observations for PSR J0108–1431 exist between 400 MHz to 1.4 GHz and the average profile has a single component across these frequencies. The outer half power width of the profile at 430 MHz is 18${}^{\circ}\pm 1^{\circ}$ (Lorimer, 1994) and that at 1370 MHz is $11\aas@@fstack{\circ}5\pm 0\aas@@fstack{\circ}3$ (Johnston & Kerr, 2018). These values are consistent with the phenomenon of radius to frequency mapping (RFM) in pulsars wherein the pulse widths are seen to decrease with increasing frequency. High-quality average polarization observations at 1.37 GHz by Johnston & Kerr (2018) reveal that the pulsar is highly polarized with about 76% linear and 15% circular polarization. The polarization position angle (PPA) shows a smooth, flatter, and almost linear traverse across the pulse profile (see also Figure 11, where the average profile from the data from Johnston & Kerr (2018) is plotted).

Based on the available data, the Empirical Theory of pulsar emission (ET:I-VII, Rankin 1983a, b, 1986, 1990; Radhakrishnan & Rankin 1990; Rankin 1993; Mitra & Rankin 2002) may be used to infer the geometry of the system. The flat PPA traverse is consistent with the rotating vector model, according to which the electric field of the radio emission lies either in the plane or perpendicular to the dipolar magnetic field line planes. As a result, the PPA traces the magnetic field planes, and near the region where the line of sight crosses the fiducial plane containing the rotation and the magnetic axis, a characteristic S-shaped curve is seen. The shape of the S-curve is a function of the angle between the rotation axis and magnetic axis $\alpha$ and the angle between the magnetic axis and the observer’s line of sight $\beta$. The slope of the PPA is maximum near the profile center ($\mid d\Psi/d\phi\mid_{max}$) and is related to the geometrical angles as,

[TABLE]

The value of the slope is estimated to be $\sim 1.2$ for J0108. In principle, the PPA traverse can be fitted to obtain the values of $\alpha$ and $\beta$, however, the fitted $\alpha$ and $\beta$ are highly correlated and they do not provide any meaningful values (von Hoensbroech & Xilouris 1997; Everett & Weisberg 2001; Mitra & Li 2004).

Instead, we can use the ET theory, according to which the pulsar radio emission beam is in the form of a core-cone structure, where, in the vicinity of the dipolar magnetic pole, there is a Gaussian-like core emission surrounded by two nested conal emission, namely the outer and inner cone. The beam is roughly circular in shape having a beam radius $\rho$, which is related to $\alpha$, $\beta$, and pulse width $\phi$ as,

[TABLE]

In ET:IV, it was found that two distinct values of $\rho$ exist at 1 GHz for the outer and inner conal emission, which are related to the pulsar period as $\rho_{outer}=5.7\,P^{-0.5}$ and $\rho_{inner}=4.3\,P^{-0.5}$. Further, in ET:VII, it was found that RFM is associated with outer conal emission, and since RFM is seen in J0108, we will assume here that the observed emission is from outer conal emission.

We calculate $\rho_{outer}=6.3^{\circ}$ using J0108’s period $P=0.807$ s. Now we can insert the value of $\rho_{outer}$, and $\beta$ obtained in terms of $\alpha$ from eqn. (4), in eqn. (5) and use an iterative procedure to find appropriate values of $\alpha$ and $\beta$ that will reproduce the observed 1 GHz half-power pulse width.Following this method we found $\alpha\sim 7.5^{\circ}$ and $\beta\sim 6.2^{\circ}$. From $\beta/\rho\sim 0.98$, we infer that the line of sight traverses the beam rather tangentially, and since the pulse profile has a single component, we can classify the pulsar to be a conal single as per the classification scheme of ET. Assuming the neutron star radius to be 10 km, a star-centered dipolar structure for its magnetic field, and same emission height across the radio profile, we can compute the radio emission height $h=10\,P\,(\rho/1.23)^{2}\sim 211$ km.

Our analysis above suggests that the radio emission properties of PSR J0108–1431 are consistent with emission arising from regions of open dipolar magnetic field lines, a few hundred km above the stellar surface. Since the line of sight cuts the emission beam rather tangentially, as suggested by the large $\beta$ values, we classify the profile as a conal single. The peak of the profile hence can be considered to lie in the fiducial plane containing the rotation and the dipolar magnetic axis.

6.2 Thermal X-ray and radio profile offset

In this section, we discuss the method to measure the offset between the thermal X-ray emission peak and the radio profile peak. We select a few X-ray events near the phase of $0.15-2$ profile peak as our X-ray times-of-arrival (TOAs). We produced radio TOAs from archival Parkes profile data made available through the CSIRO data access portal. Obtaining a combined timing solution, while allowing a constant offset parameter between the radio and X-ray TOAs, allows us to determine the absolute profile peak offsets between the two energy bands.

We perform non-parametric bootstrap re-sampling (Feigelson & Babu, 2012) of the original list of event phases to obtain 10000 samples of event phases in the $0.15-2$ keV range. The 10000 KDE-smoothed profiles produced from these sampled event phases are expected to represent the same underlying distribution (the true pulse profile), varying only due to the statistical variance inherent in the data (observed profile). The peak of the $0.15-2$ keV profile is estimated to be at $\tilde{\phi}_{x}=0.003(40)$, where, the value in the parenthesis represent the uncertainty at the $10-90$ percentile level from the quoted median value. We use the barycentered arrival times in MJD of eight randomly selected events from the narrow $0.003\pm 0.003$ phase range to produce the list of X-ray TOAs, assigning a 33 ms uncertainty (corresponding to $\sigma_{\phi}=0.04$) to each TOA. Due to low S/N, neither the time sub-integrated data produce reliable profiles nor the integrated profile produces a good template. Hence, the X-ray events near the original profile peak and the peak uncertainty estimated through bootstrap re-sampling together produce the most conservative TOAs for the given data.

Archival data from the Parkes telescope for J0108 is available as part of a regular monitoring campaign of pulsars with high energy emission (Weltevrede et al., 2010). We obtained radio TOAs over a long baseline (April 2009 to September 2011 $\equiv$ MJD $54931-55821$) around the X-ray observation epoch (June 2011 $\equiv$ MJD 55728). We used PSRCHIVE (Hotan et al. 2004; van Straten et al. 2012) to clean and reduce the sub-banded and sub-integrated profile data and produce integrated profiles. We created a noise-free template from a very high S/N profile, and finally, obtain the TOAs using the Fourier phase gradient algorithm (Taylor, 1992).

We fit the TOAs using TEMPO2 (Hobbs et al., 2006) through the python wrapper libstempo555http://vallis.github.io/libstempo/. Extrapolating the pulsar’s timing solution to the epoch of X-ray observations require updates to only the rotation parameters: frequency ($F0$), frequency derivative ($F1$), and the frequency second derivative ($F2$). We include the X-ray TOAs and introduce a constant offset parameter (JUMP) to find the radio-X-ray offset (Table 5; Figures 12 and 13). The system delays in the Parkes radio data are of the order of a micro-second (Manchester et al., 2013) which is orders of magnitude lower than our TOA uncertainties and hence are ignored. The offset of the X-ray peak from the radio peak is estimated to be $\Delta\phi_{\rm x-r}=-0.04947(36)$ where the value in the parenthesis represent the uncertainty at the $10-90$ percentile level from the quoted median value. Subtracting this from the phase of X-ray peak, the phase of the radio peak is estimated to be $\phi_{\rm r}=0.037_{-0.059}^{+0.041}$.

In order to obtain the peak of the thermal emission, we fit a sinusoidal model curve to the $0.15-0.7$ keV profile, in the $0.2-0.7$ phase range where the emission is thermal (Figure 14). The model curve is of the form

[TABLE]

where, $A0,A$, and $\phi_{0}$ are the constant level, the amplitude of the Sine, and the zero phase, respectively. We obtain a Bayesian 3 parameter fit, using MCMC to sample the posterior with 100 walkers taking 3000 steps after a burn-in of 500 steps. Our prior for the zero-phase of the sine curve ($\phi_{0}$) is informed by the constraint on the thermal peak obtained from spectral fitting in Section 5.1. The posterior distribution of the peak phase has 98% of the values spread in the $0.2-0.7$ phase range and only 2% outside. Hence, we assign a probability of 0.98 that Sine curves peak in the $0.2-0.7$ phase range and 0.02 elsewhere. The rest of the parameters are assigned uniform priors. The peak phase of thermal emission is estimated to be $\tilde{\phi}_{\rm th}=0.43\pm 0.14$ where, the uncertainties are at the $10-90$ percentile levels. From the distribution for peak of the radio and thermal X-ray emission, we estimate the former to lead the latter by $\Delta\phi_{\rm r-th}=0.41_{-0.20}^{+0.13}$ (Figure 15). We also estimate a 99.7% probability that the offset between the thermal X-ray peak and the radio peak is greater than 0.1. This corresponds to a shift of $S\sim 2\,\pi\,\Delta\phi\,R_{\rm NS}\,\sin\alpha\gtrsim 820$ between the hotspot and the dipole axis.

7 Discussions and Conclusions

The need for a multipolar surface magnetic field has been recognized early on to explain significant pair multiplication above the pulsar polar cap (Ruderman & Sutherland, 1975). The presence of a strong multipolar magnetic field is often invoked to explain X-ray thermal emission from areas that are apparently smaller than the polar cap area for a dipolar surface magnetic field. Recent compilations of polar cap radii from X-ray observed non-recycled pulsars have been used as strong observational evidence for a multipolar surface magnetic field (Geppert 2017; Rigoselli & Mereghetti 2018). We revisit the literature listed in these compilations and assess the robustness of their thermal polar cap area measurements.

We first identify systems with definite evidence of thermal emission. Of the 18 pulsars in the lists, 9 of the pulsars do not require a thermal model component to explain the spectra (PSR J0108$-$1431; Posselt et al. 2012, PSR B1133+16; Szary et al. 2017, PSR B1451$-$68; Pancrazi et al. 2012, PSR B0628$-$28; Tepedelenlıoǧlu & Ögelman 2005, PSR B0834+06; Gil et al. 2008, PSR B1719$-$37; Oosterbroek et al. 2004, PSR B1929+10; Misanovic et al. 2008, PSR J2043+2740; Becker et al. 2004, PSR B2224+65; Hui et al. 2012). For these pulsars, either the S/N of the data is too low to differentiate thermal contribution from the apparent single-component non-thermal emission, or a single PL component adequately explains the spectra. In such cases, thermal polar cap areas are often estimated by over fitting the spectra with extra BB components, fixing the ISM absorption or non-thermal model parameters to obtain BB parameters, or from an upper-limit on the thermal component flux. The PC area estimates from such systems are unreliable for scientific inference.

The estimates of the thermal emission area found in literature depends crucially on the assumption that the emission is a BB. So, among the remaining 9 pulsars, none of them show definitive evidence for isotropic emission of a BB. On the contrary, for pulsars with high S/N data, the pulse profiles are sharper and non-sinusoidal in the energy range where the polar cap emission dominates often with high pulsed fractions or increasing pulsed fraction with energy (PSR B0656+14; De Luca et al. 2005; Arumugasamy et al. 2018, PSR B1055$-$52; De Luca et al. 2005, PSR J0633+1746; Kargaltsev et al. 2005). These are indicative of beamed (anisotropic) emission significantly deviating from a BB. For rest of the pulsars, either no tests were performed and/or no evidence found to prove the existence of a BB (B0943+10; Mereghetti et al. 2016, PSR B1822$-$09; Hermsen et al. 2017, PSR B0355+54; Klingler et al. 2016, PSR B0114+58; Rigoselli & Mereghetti 2018, PSR J1740+1000; Kargaltsev et al. 2012, PSR J2021+4026 Lin et al. 2013).

Once the assumption of the BB emission model is relaxed, recovering the polar cap area from the observed thermal emission is not straight-forward. In fact, as shown in this paper, fitting a BB to beamed emission from a neutron star atmosphere leads to under-estimate of the polar cap area. BB area estimates found in literature also do not fully propagate the distance uncertainties which for DM based distances are usually not well quantified and for some parallax based measurements can be significant. Hence, the observed PC radius estimates under the BB assumption are reliable only if there is clear evidence of thermal emission from the NS surface, followed by evidence that the thermal emission is a BB and not processed and beamed, and finally, the uncertainties in the distance and other component parameters are included.

We evaluate the evidence of a global dipolar magnetic field structure for the very old pulsar PSR J0108–1431 in this work. Our results emphasize how low S/N data is incapable of uniquely determining the thermal emission model. We answer the following questions:

1. (A)

Does the phase-integrated spectrum show clear evidence of a thermal component?

A single $\Gamma\approx 3$ PL component is sufficient to explain the phase-integrated spectrum. A single BB is not a good fit and two-component models (PL+BB or PL+NSA) over-fit the spectra, and hence are not required. We use the results of a PL+BB fit to get the rough energy ranges in which either of the two components could dominate. However, parameter estimates from such over-fit two-component models are unsuitable for any objective scientific inference. 2. (B)

Is there evidence of multiple emission components (including thermal) in pulse profiles or phase-separated spectra?

We see a strong indication of two distinct emission components separated in phase when comparing the profiles extracted in the $0.15-0.7$ keV and $0.7-2$ keV ranges. The softer component, seen in the phase range $0.2-0.7$, is only fit with a thermal component (BB or NSA) or an unrealistically steep, $\Gamma\approx 5$ PL. The emission in the rest of the phase range ($0.7-0.2$) is best fit with a PL (BB or NSA models fit poorly) and hence predominantly non-thermal. 3. (C)

Is the soft emission component fit only by a blackbody model?

The soft spectrum in the $0.2-0.7$ phase range excludes a non-thermal description via a PL model. The spectrum can be explained with a 120 eV BB or a 78 eV NSA model emitting from a fraction, $3\times 10^{-5}$ or $3\times 10^{-4}$ of the NS surface, respectively. The statistics are insufficient to distinguish between the two models. 4. (D)

Can the PC area be measured reliably?

For an observer at infinity (far from the gravitational influence of the NS), the polar cap area for a star-centered, spin-aligned, global dipolar magnetic field rotating with a period of $0.808\,$s will appear to be $A^{\infty}_{\rm d,pc}=0.139$ km2 (assuming $M_{\rm NS}=1.4\,M_{\sun}$ and $R_{\rm NS}=10$ km). The median BB emitting area estimate is a factor of 2 lower, with 86% of the distribution below the dipolar value, after including the distance uncertainties. However, the median NSA emitting area estimate is a factor of 4 higher, with 80% of the distribution above the dipolar value. Since the two models are statistically indistinguishable, we conclude that the area of the polar cap cannot be measured reliably with the current models. 5. (E)

Are there any additional results from the X-ray analysis?

We found an absorption-like feature in the spectrum extracted from half of the rotational phase (thermal continuum in the $0.2-0.7$ phase range). The feature cannot be explained by instrumental systematics or ISM absorption because it is absent in the spectrum from the rest of the phase. We fit the spectral feature with a Gaussian optical depth profile to obtain central absorption energy $E_{\rm c}=0.33_{-0.08}^{+0.05}$ keV with equivalent width $EW=0.105\pm 0.05$. The absorption energy corresponds to local magnetic field strengths of $(2.2-3.3)\times 10^{10}$ G or $(4-6)\times 10^{13}$ G assuming electron or proton cyclotron absorption, respectively. This could be evidence for proton cyclotron absorption in the NS atmosphere with surface magnetic field strength exceeding $\sim 10^{13}$ G. Such features with multiple harmonics and/or phase-dependent variation have been confirmed in high magnetic field pulsars (e.g., Sanwal et al. 2002; De Luca et al. 2004; Tiengo et al. 2013; Borghese et al. 2015) where they are attributed to cyclotron absorption by protons confined in multipolar magnetic field structures close to the NS surface. 6. (F)

What are the constraints on the coherent radio emission region?

In section  6 we have demonstrated that the coherent radio emission from J0108 is consistent with the model of emission arising from a circular beam with underlying open dipolar magnetic field line region, at about few hundred km above the neutron star surface. We also infer that the pulsar is an almost aligned rotator, with the observer line-of-sight cutting the emission beam tangentially. For a star centered static magnetic dipole, the thermal hot spot and the radio emission region should be aligned in time. However, for the rotating case at an emission height of $\sim 211$ km, aberration and retardation effects (Blaskiewicz et al. 1991; Dyks et al. 2004) advance the radio emission phase by $\Delta\phi\sim 0.011$, while magnetic sweepback (Dyks & Harding, 2004) introduces a lag $\Delta\phi\sim 0.015$ with respect to the phase of thermal hotspot emission. Hence, the radio peak is expected to lag the X-ray peak by $\Delta\Phi\sim 0.004$ for a rotating dipolar magnetic field. 7. (G)

Is the inferred surface thermal peak and radio peak offset for J0108 consistent with a star-centered dipole?

We estimate an offset $\Delta\phi\approx 0.4$ between the X-ray thermal peak and the leading radio peak, with a 99.7% probability that the offset is greater than 0.1. Thus, the dislocation of the thermal polar cap with respect to the magnetic dipolar axis is $S\sim 2\,\pi\,\Delta\phi\,R_{\rm NS}\,\sin\alpha\gtrsim 0.8$ km, assuming $R_{\rm NS}=10$ km (see also Szary et al. 2017). An offset of such magnitude cannot be explained by a star-centered global dipolar magnetic field configuration. Hence, in the case of J0108, the measured X-ray thermal and radio emission offset provides a strong evidence for multipolar surface magnetic fields.

We choose PSR J0108–1431 as a prototypical system to show how the prevailing method of determining PC area to infer the presence of a multipolar surface magnetic field is unreliable and how the radio and X-ray thermal emission offset method provides a better alternative. In follow-up work, we plan to extend such analyses to the remaining sample of radio RPPs with thermal PC emission and determine whether strong evidence for multipolar magnetic fields can be found with the currently available data.

Acknowledgements

We appreciate the ATNF and CSIRO Data Access Portal for making such high-quality data publicly available. We thank the XMM-Newton help-desk for important clarifications regarding the XMM-Newton data analysis. We also thank Bettina Posselt, Oleg Kargaltsev, George G. Pavlov, and Joanna Rankin for the useful discussions. We greatly appreciate Dr. Sebastien Guillot for carefully reviewing our manuscript and suggesting essential improvements. DM acknowledges funding from the grant “Indo-French Centre for the Promotion of Advanced Research - CEFIPRA”.